A0-A113  998  AZR  FORCE  INST  OF  TECH  HR16HT-PATTERS0N  AFB  OH  SCHOO— ETC  P/6  12/1 

A  UNIVARIATE  MONTE  CARLO  TECHNIQUE  TO  APPROXIMATE  RELIABILITY  C->ETC(U) 
dec  79  R  B  PUTZ 

UNCLASSIFIED  AFIT/S0R/MA/79D-7  NL 


UNITED  STATES  AIR  FORCE 
AIR  UNIVERSITY 

AIR  FORCE  INSTITUTE  OF  TECHNOLO 

Wright<*Pott«rton  Air  Forc«  Bat«,Ohio 


ELECT 


AFIT/GOR/MA/79D-7 


d) 


A  UNIVARIATE  MONTE  CARLO  TECHNIQUE  TO 
APPROXIMATE  RELIABILITY  CONFIDENCE 
LIMITS  OF  SYSTEMS  WITH  COMPONENTS 
CHARACTERIZED  BY  THE 
WEIBULL  DISTRIBUTION 

THESIS 


AFIT/GOR/MA/79D-7  Randall  B.  Putz 

Captain  USAF 


Approved 


for  public 


release ; 


distribution  unlimited 


t 


AFIT/GOR/MA/79D-7 


A  UNIVARIATE  MONTE  CARLO  TECHNIQUE  TO  APPROXIMATE 
RELIABILITY  CONFIDENCE  LIMITS  OF  SYSTEMS  WITH 
COMPONENTS  CHARACTERIZED  BY  THE 
WEIBULL  DISTRIBUTION 

THESIS 


Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 
Air  University 

in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 
Master  of  Science 


by 


Randall  B.  Futz,  B.S. 
Captain  USAF 


December  1979 


Approved  for  public  release;  distribution  unlimited 


Preface 


Reliability  engineering  is  a  critical  area  in  the 
design  of  any  system.  Too  often  a  promising  and  exciting 
system  which  is  on  the  threshold  of  technology  turns  into 
a  nightmare  to  maintain  because  it  is  constantly  breaking 
down.  Improvements  are  clearly  needed  in  the  field  of 
estimating  system  reliability.  For  several  years,  thesis 
students  under  the  direction  of  Dr.  Albert  H.  Moore  have 
been  doing  research  in  the  area  of  reliability,  especially 
with  regard  to  Monte  Carlo  simulation  of  systems  reliabil¬ 
ity.  This  thesis  explores  the  possibility  of  simplifying 
reliability  analysis  while  still  maintaining  reasonable 
accuracy . 

I  want  to  thank  Dr,  Moore  for  his  patient  instruc¬ 
tion  over  the  last  months  and  his  continuing  Interest  in 
this  thesis.  I  also  would  like  to  thank  my  wife  Nancy  for 
the  daily  encouragement  she  has  given  me  for  the  past 
eighteen  months. 

I  would  like  to  thank  Phyllis  Reynolds  for  typing 
my  thesis. 
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Abstract 


A  univariate  Monte  Carlo  technique  Is  developed  for 
the  determination  of  lower  confidence  limits  of  system  reli¬ 
ability  based  on  component  test  data.  It  Is  assumed  that 
the  component  teat  data  consists  of  failure  times  which  are 
distributed  according  to  a  known  two-parameter  Weibull  prob¬ 
ability  distribution.  These  failure  times  are  randomly 
generated  using  the  true  shape  and  scale  parameters  of  the 
distribution.  Maximum-likelihood  estimators  are  found  for 
the  shape  and  scale  parameters  and  then  substituted  Into 
the  reliability  equation  to  obtain  the  maximum- likelihood 
estimator  for  the  component  reliability.  The  estimated  bias 
in  this  estimator  is  subtracted  to  yield  an  approximately 
unbiased  estimator  of  the  component  reliability.  Using  the 
empirical  variance  of  the  reliability  estimate  and  assuming 
a  normal  distribution,  a  Monte  Carlo  simulation  is  run  for 
four  hypothetical  systems  consisting  of  as  many  as  five  com¬ 
ponents.  The  simulation  is  repeated  600  times.  Since  the 
true  reliability  is  known,  on  each  run  it  can  be  determined 
if  the  desired  confidence  intervals  contain  the  true  system 
reliability.  The  result  is  an  absolute  measure  of  the 
effectiveness  of  the  univariate  technique.  The  entire  simu¬ 
lation  was  run  for  component  test  data  sample  sizes  ranging 
from  ten  to  one  hundred.  A  second  run  of  600  was  made  to 
examine  the  Monte  Carlo  sampling  error  for  component  test 
sample  sizes  of  100. 
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A  UNIVARIATE  MONTE  CARLO  TECHNIQUE  TO  APPROXIMATE 


RELIABILITY  CONFIDENCE  LIMITS  OF  SYSTEMS  WITH 
COMPONENTS  CHARACTERIZED  BY  THE 
WEIBULL  DISTRIBUTION 

I .  Introduction 

Problem  Statement 

Testing  a  complex  system  In  order  to  determine  Its 
reliability  can  be  expensive  In  terms  of  time  and  money. 

Even  the  preliminary  design  of  a  system  requires  advance 
knowledge  about  how  a  network  of  components  will  have  to 
be  structured  to  yield  the  desired  system  reliability.  The 
purpose  of  this  thesis  is  to  determine  system  reliability 
based  on  reliability  data  about  individual  components 
within  a  system. 

Reliability 

The  reliability  of  a  system  is  defined  as  the  proba¬ 
bility  that  the  system  will  perform  its  intended  function 
for  a  specified  period  of  time  under  stated  conditions. 

The  current  importance  of  reliability  is  due  in  large  part 
to  the  experience  of  the  U.S.  military  during  and  immedi¬ 
ately  after  World  War  II.  Tremendous  difficulties  were 
encountered  in  assuring  reliable  equipment  for  the  fighting 
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forces.  For  example,  during  Che  war,  50%  of  all  stored  air¬ 
borne  electronic  equipment  became  unserviceable  before  even 
being  used.  The  Army  reported  high  truck  engine  and  power 
plant  mortality  rates.  In  1947,  the  Navy  stated  that  70% 
of  its  electronic  equipment  did  not  operate  satisfactorily 
(Ref  12:Sec  I,  1-2).  Although  reliability  theory  and  tech¬ 
niques  for  reliability  estimation  have  shown  magnitudes  of 
improvement  since  World  War  II,  there  is  still  the  need  for 
significant  developments.  As  systems  become  more  complex 
and  testing  becomes  more  expensive  (look  at  the  space 
shuttle  as  an  example) ,  accurate  reliability  estimation  tech¬ 
niques  based  on  small  sample  sizes  must  be  developed. 

The  Welbull  Distribution 

The  Weibull  probability  density  function  was  ini¬ 
tially  developed  for  the  study  of  fatigue  failure  in 
materials.  Later,  it  was  discovered  that  electron  tubes  of 
various  types  fail  according  to  the  distribution  (Ref  12: 

Sec  II,  55).  Thus  the  Weibull  can  be  used  to  express  the 
distribution  of  time  to  failure  for  many  mechanical  and 
electrical  devices.  In  addition,  there  is  frequently  the 
assumption  that  component  failures  (especially  electronic 
components)  are  distributed  according  to  the  exponential 
distribution.  Zelen  and  Dannemlller  (Ref  15)  have  pointed 
out  that  the  exponential  distribution  is  generally  not  a 
robust  approximation  to  the  Weibull,  especially  if  the  shape 
parameter  is  greater  than  1. 
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The  Welbull  density  function  Is 


f(t;k,e,c)  =  ■  ^  e  k,9  ^  0;  c  <  t  (1) 

e*' 

=■  0  elsewhere, 

where  t  Is  the  time,  6  Is  the  scale  parameter,  c  Is  the 
location  parameter,  and  k  Is  the  shape  parameter.  The 
scale  parameter  determines  the  spread  or  dispersion  of  the 
function  about  Its  mean.  The  location  parameter  determines 
the  point  of  origin  of  the  distribution.  The  shape  param¬ 
eter  determines  whether  the  failure  rate  is  decreasing, 
increasing,  or  Invariant  with  time.  This  gives  the  Welbull 
distribution  flexibility  to  model  components  with  failure 
rates  that  decrease,  increase,  or  remain  constant  over  time. 
For  example,  if  the  shape  parameter  is  1,  the  Weibull 
degenerates  to  the  exponent ial  density  function.  If  the 
shape  parameter  is  3.5,  the  Weibull  is  an  excellent  approxi¬ 
mation  to  the  normal  distribution.  A  shape  parameter  of 
2  yields  a  function  which,  when  correctly  scaled,  can  approx¬ 
imate  Beta  distributions  that  are  skewed  towards  the  right. 
Fig.  1  shows  the  flexibility  of  the  Weibull  to  model  the 
exponential,  the  normal,  and  the  Beta  distributions. 

An  individual  component  that  follows  a  Weibull  dis¬ 
tribution  in  time  to  failure  will  have  its  reliability 
expressed  by 

R(t)  =  exp(-((t-c)  /e)'")  .  (2) 
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Fig.  1.  Weibull  Distribution  with  Shape  Parameters 

of  1 ,  2 ,  and  3 . 5 

Therefore,  determining  the  reliability  of  one  component  is 
a  relatively  easy  procedure  of  substituting  into  the  reli¬ 
ability  equation.  However,  one  is  frequently  interested 
in  determining  the  reliability  of  a  system  of  several  com¬ 
ponents,  each  possessing  a  time  to  failure  given  by  the 
Weibull  probability  density  function.  Problems  arise  in 
that  analytical  techniques  designed  to  solve  for  the  reli¬ 
ability  of  complex  systems  quickly  become  overwhelming  as 


the  number  of  components  increases.  Therefore,  computer 


simulation  has  become  a  handy  tool  to  utilize  on  complex 
systems . 

Survey  of  Past  Work 

Some  of  the  original  work  in  the  use  of  Monte  Carlo 
methods  for  determining  confidence  limits  on  system  relia¬ 
bility  was  done  by  Donald  Orkand  (Ref  8)  in  1960.  He 
pointed  out  the  serious  limitations  in  using  only  point 
estimates  of  component  reliability  to  determine  a  point 
estimate  of  the  system  reliability.  He  also  showed  the 
flexibility  of  using  the  Monte  Carlo  method  and  its  adapta¬ 
bility  to  the  circuitry  of  various  complex  systems.  Albert 
Bernhoff  (Ref  2)  explained  the  problems  of  applying  ana¬ 
lytical  approaches  to  find  system  reliability  confidence 
intervals  by  combining  component  reliability  confidence 
intervals.  He  concluded  that  the  analytical  approaches 
generally  do  not  work.  In  1967,  Louis  Levy  and  Albert  Moore 
(Ref  5)  developed  a  Monte  Carlo  technique  for  obtaining 
reliability  confidence  intervals  for  systems  with  components 
characterized  by  the  exponential,  normal,  lognormal.  Gamma, 
and  Weibull  distributions.  Robert  Lannon  (Ref  4)  estab¬ 
lished  reliability  confidence  intervals  for  the  Weibull  dis¬ 
tribution  using  a  bivariate  analysis  of  the  shape  and  scale 
parameters.  Robert  Snead  (Ref  11)  studied  reliability 
estimates  using  the  Weibull,  Gamma,  and  Logistic  distribu¬ 
tions.  For  each  distribution,  Snead  used  the  property 
that  the  reliability  estimates  are  asymptotically  normal. 
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Thus  he  was  able  Co  perform  a  univariate  analysis  using 
Just  Che  reliability  parameter,  as  long  as  large  sample 


sizes  existed.  However,  because  Che  primary  value  In  the 
Monte  Carlo  method  Is  to  reduce  the  prohibitively  expensive 
testing  of  large  samples  of  systems  and  their  individual 
components,  it  would  be  extremely  helpful  if  a  univariate 
method  could  be  developed  for  small  sample  sizes. 

Darrel  Thoman,  Lee  Bain,  and  Charles  Antle  (Refs 
13  and  14)  have  done  considerable  research  Into  Che  Welbull 
probability  distribution  and  the  maximum  likelihood  esti¬ 
mation  of  the  Welbull  parameters.  Key  results  from  their 
studies  show  that,  assuming  the  two-parameter  Welbull  (that 
is,  the  shape  and  scale  parameters  are  unknown  while  the 
location  parameter  is  zero),  the  distribution  of  the  maxi- 
mum  likelihood  estimator  of  reliability,  R(t),  depends  only 
on  the  true  system  reliability  R(t)  and  the  sample  size  n. 
Also,  R(t)  is  nearly  unbiased  and  has  a  variance  practi¬ 
cally  equal  to  the  Cramer-Rao  lower  bound  for  the  variance 
of  an  unbiased  estimator  (Ref  13:363).  These  results  form 
the  basis  for  the  research  done  in  this  thesis. 

Ob.1  ec  t  ives 

There  are  two  objectives  of  this  study.  The  first 
is  to  simplify  the  bivariate  analysis  of  Lannon  by  mapping 
the  shape  and  scale  parameters  onto  the  reliability  param¬ 
eter,  resulting  in  a  univariate  analysis.  The  second  is  Co 
explore  how  far  one  can  reduce  the  sample  sizes  ot  component 
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test  data  while  retaining  reasonably  accurate  confidence 
Intervals . 

Approach 

Several  assumptions  are  necessary  In  order  to  make 
the  problem  tractable.  The  first  is  that  system  components 
have  failure  times  represented  by  the  Weibull  probability 
distribution.  Also,  all  components  must  have  been  sub¬ 
jected  to  reliability  testa,  allowing  estimation  of  the 
shape,  scale,  and  location  parameters.  The  distribution  of 
the  estimators  for  the  parameters  must  also  be  known. 
Finally,  there  must  be  a  known  mathematical  relationship 
between  the  system  reliability  and  component  reliabilities 
(Ref  7:459).  In  this  thesis,  it  is  assumed  that  the  com¬ 
ponents  fail  independently  of  one  another.  This  allows  for 
the  use  of  standard  formulas  for  the  reliability  of  systems 
having  series  or  parallel  connections,  and  for  the  reduction 
of  complex  systems  into  combinations  of  series  and  parallel 
configurations.  Orkand  (Ref  8:7-8)  describes  procedures 
to  apply  if  there  is  a  dependence  relation  among  system 
component  s . 

The  fundamental  approach  that  will  be  applied 
throughout  this  thesis  is  to  determine  reliability  esti¬ 
mates  and  the  associated  confidence  intervals  using  the 
Monte  Carlo  method.  This  method  uses  random  sampling  to 
investigate  the  solution  of  deterministic  or  stochastic 
problems.  In  essence,  one  inputs  a  set  of  random  variables 
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from  specified  probability  distributions,  and  then  solves 
the  problem  for  each  set  of  Inputs  to  obtain  a  random 
sample  of  outcomes  (Ref  7:459).  The  steps  to  the  method 
are  as  follows: 

1.  Using  the  component  test  data  (set  of  failure 
times),  calculate  the  maximum  likelihood  estimators  for  the 
parameters  of  the  underlying  distribution.  Do  this  for 
all  components. 

2.  Determine  the  Joint  distribution  of  the  maximum 
likelihood  estimators.  Repeat  this  for  all  components. 

3.  Generate  a  sample  value  for  each  parameter  by 
sampling  from  the  distribution  of  the  parameter. 

4.  Substitute  the  sample  parameters  into  the  com¬ 
ponent's  reliability  equation  to  derive  a  sample  component 
reliability . 

5.  Repeat  steps  3  and  4  for  each  component,  and 
then  substitute  the  sample  component  reliabilities  into  the 
system  reliability  equation  to  obtain  a  sample  system  reli¬ 
ability  . 

6.  Repeat  steps  3-5  to  obtain  the  desired  number 
of  samples  of  the  system  reliability. 

7.  Order  the  system  reliability  samples  and  deter¬ 
mine  the  desired  confidence  limits. 

The  primary  advantage  of  the  Monte  Carlo  method  is 
its  adaptability  to  any  type  of  complex  system.  It  is  a 
proven  method  that  can  be  used  for  the  initial  investigation 
into  any  systems  reliability  problem  that  fits  the 
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above-mentioned  assumptions.  It  Is  also  a  valuable  method 
to  apply  when  other  methods  fall  due  to  the  complexity  of 
the  system  (Ref  7:459). 
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II .  Theoretical  Development 

Step  2  of  the  Monte  Carlo  method  requires  knowledge 
about  the  joint  distribution  of  the  maximum  likelihood 
estimators.  Since  an  objective  of  this  thesis  is  to 
research  a  univariate  Monte  Carlo  method,  the  distribution 

/V 

of  the  single  estimator,  R(t),  is  Important. 

Maximum  Likelihood  Estimation 

An  estimator  is  simply  an  approximation  to  the 
unknown  value  of  some  population  parameter  (Ref  9:39).  For 
example,  if  one  is  interested  in  the  mean  utility  bill  in 
some  city,  he  can  sample  from  several  different  households 
and  use  the  sample  mean  as  an  estimator  to  the  mean  utility 
bill  for  all  the  households  within  the  city.  Also,  the 
sample  standard  deviation  is  frequently  used  as  an  esti¬ 
mator  for  the  population  standard  deviation.  A  particular 
type  of  estimator  is  the  maximum  likelihood  estimator  (MLE)  . 
An  example  can  best  explain  this  estimator.  Assume  that 
someone  has  a  box  containing  ten  balls,  each  of  which  is 
either  black  or  white,  and  is  asked  to  estimate  the  number 
of  black  balls  based  on  a  sample  of  three  balls  drawn  from 
the  box.  If  the  sample  is  three  black  balls,  then  the 
maximum  likelihood  estimator  of  the  total  number  of  black 
balls  is  ten.  This  is  because  a  population  of  ten  black 
balls  would  maximize  the  probability  of  drawing  the  sample 
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of  three  black  balls.  In  this  case,  that  probability  would 

be  one.  Maximum  likelihood  estimators  are  those  values  of 

the  parameters  which  maximize  the  probability  or  the  joint 

density  (the  likelihood)  of  the  observed  sample  (Ref  6:302). 

Some  symbols  and  special  terminology  need  to  be 

introduced  here.  Let  *^2  ' '  *  * ’^n  sample  observations 

taken  on  the  corresponding  random  variables  Y,,Y„ . Y  . 

1  z  n 

Then  if  Yt,Y_,...,Y  are  discrete  random  variables,  the 
j.  z  n 

likelihood  of  the  sample,  L,  is  defined  to  be  the  joint 

probability  of  yi,yo,...,y  .  If  Y.,Y.,,..,,Y  are  continuous 

1  z  n  j.  z  n 

random  variables,  then  the  likelihood  L  is  defined  to  be 

the  joint  density  evaluated  at  yj^»y2 . y^j*  either  the 

discrete  or  the  continuous  case,  maximum  likelihood  esti¬ 
mators  are  those  values  of  the  estimators  which  maximize  L, 
the  likelihood  of  the  sample  (Ref  6:303). 

In  researching  the  reliability  of  systems  with  com¬ 
ponents  that  fail  according  to  the  Weibull  probability  dis¬ 
tribution,  the  sample  is  a  set  of  n  failure  times  (say  t^, 
t2,...,t^)  obtained  by  testing  each  component.  The  objec¬ 
tive  is  to  pick  the  values  of  the  shape  and  scale  parameters 
(location  parameter  is  assumed  to  be  zero)  that  will  maxi¬ 
mize  the  likelihood  of  the  observed  failure  times.  As 
defined  earlier,  the  likelihood  L  equals  the  joint  density 
of  the  sample. 

L  »  f  (tj^,t2».  .  •  k,0) 
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But  since  the  failure  times  are  Independent  of  one  another 


L  -  f(t,;  k,0)f(t„;  k,0)...f(t  ;  k,0) 
1  ^  n 


To  maximize  L,  the  usual  procedure  Is  to  take  the  partial 
derivatives  of  L  with  respect  to  the  parameters,  set  the 
derivatives  equal  to  zero,  and  then  solve  for  the  values  of 
the  parameters.  But  maximizing  ln(L)  Is  easier  than  maxi¬ 
mizing  L  and  will  result  In  the  same  estimators  for  the 
shape  and  scale  parameters. 

n 

ln(L)  -  ln[  TT  f(t  ;  k.0)] 

1-1  ^ 


In 


1-1 


kt 


k-1 


kl 


n 

-  ^[ln(k)  +  (k-l)ln(t^)  -  k  In  (0) 
1-1 


-  n  In  k  + 


(k-1) 


In  t 


1-1 


nk  In  0 


(3) 


Now  taking  the  partial  derivative  of  ln(L)  with  respect 
to  0, 
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3  In  L 
39 


-nk/0 


+  ke 


-k-1 


i-1 


(4) 


Setting  Eq  (4)  equal  to  zero  and  solving  for  6  yields 


/N 
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(5) 


Substituting  the  value  of  0  Into  Eq  (3)  and  taking  the  par¬ 
tial  derivative  of  ln(L)  with  respect  to  k. 


3  In  L 
3  k 


n/k  + 


n 


^InCti) 

1-1 


1-1  1-1 


(6) 


Setting  the  derivative  equal  to  zero  yields 


n 


1-1 


(7) 


Eqs  (5)  and  (7)  are  two  equations  In  the  two  unknowns,  k 
and  9.  Simultaneous  solution  of  these  equations  yields  the 
maximum  likelihood  estimators  of  k  and  0. 


The  Distribution  of  R(_t) 

The  MLE  of  reliability  Is  found  by  substituting  the 
MLEs  for  k  and  9  Into  the  reliability  formula,  Eq  (2). 


'  V 
'  > 


I 


I 
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(8) 


R(t)  *  e  9 

If  one  repeatedly  draws  samples  of  size  n  from  a  population 
distributed  according  to  the  Weibull,  and  uses  the  pre¬ 
viously  mentioned  technique  to  find  the  MLE  of  reliability 

A 

for  each  sample,  then  a  distribution  of  the  MLE  R(t)  will 

result.  Let  9  =»  (9/9)  and  k  *  k/k.  Thoman,  Bain,  and 

s  s 

Antle  (Ref  13)  have  proven  that  the  joint  distribution  of 

9  and  k  is  Independent  of  9  and  k.  They  go  on  to  prove 
s  s 

(Ref  14)  that  the  distribution  of  R(t)  depends  only  on  R(t). 

A. 

In  R(t)  *  -(t/9)^ 
ln[-ln  R(t)  ]  -  k  ln(t/9) 

=  k/k  ln{  (t/0  )  ^] 

In  [-ln(R  (t)  )  ]  *  tc/k  ln[  (t/9) (9/9) 

>  kg  ln[(-ln  R(t))0g~^] 

Since  k^  and  9^  are  distributed  independently  of  k  and  9, 
the  distribution  of  R(t)  depends  only  on  R(t).  Thus  the 
parameters  t,  k,  and  9  can  affect  the  distribution  or  R(t) 
only  through  R(t)  (Ref  14:363).  This  is  a  significant 
result.  It  allows  for  the  univariate  analysis  of  R(t), 
in  place  of  the  bivariate  analysis  (using  k  and  0)  that 
Lannon  did.  What  is  being  said  is  that,  if  two  components 
have  different  underlying  Weibull  distributions  (different 
shape  and  scale  parameters) ,  but  both  have  the  same 
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reliability,  then  they  will  have  Identical  distributions 
of  their  maximum  likelihood  estimators  of  reliability. 

An  important  criteria  of  an  estimator  is  its  bias. 
Usually  a  small  or  zero  bias  is  desired.  The  bias,  B,  is 
defined  as  the  expected  value  of  the  estimator  minus  the 
true  value  of  the  parameter  being  estimated.  In  this  case, 

B  -  E[R(t)]  -  R(t) 

A 

The  bias  of  R(t)  can  be  estimated  by  Monte  Carlo  simulation 
Thoman,  Bain,  and  Antle  (Ref  14:365)  derived  10,000  estl- 
r.ates  of  R(t)  using  sample  sizes  of  8  to  100  and  true  reli¬ 
abilities  of  0.5  to  0.98.  The  10,000  estimates  were  aver- 

A 

aged  to  obtain  EfRCt)].  and  then  the  true  reliability  was 
subtracted  to  yield  the  bias.  Unfortunately,  Thoman,  Bain, 
and  Antle  did  not  include  enough  significant  digits  to  be 
helpful.  David  Antoon  (Ref  1)  did  a  Monte  Carlo  analysis 
and  derived  2000  estimates  of  R(t).  The  bias  was  calcu¬ 
lated  for  the  same  range  of  sample  sizes  and  true  reliabili 
ties  as  Thoman,  Bain,  and  Antle.  Antoon's  results  compared 
favorably  and  are  shown  in  Table  I.  As  can  be  seen,  the 

A  A 

bias  of  R(t)  is  very  close  to  zero.  Thus  the  MLE ,  R(t), 
is  almost  an  unbiased  estimator  of  R(t). 

Another  desired  characteristic  of  an  estimator  is 
that  it  have  a  small  variance.  In  this  case, 

Var[R(t)]  -  ^ - 

i»l 
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TABLE 


▼ 


where  M  is  the  Monte  Carlo  size,  or  the  number  of  samples 
generated.  As  mentioned,  Thoman,  Bain,  and  Antle  generated 
10,000  samples  from  which  they  calculated  estimates  of  the 
variance  of  R(t).  Table  II  (Ref  14:366)  shows  their 
results.  As  can  be  seen,  only  three  significant  digits 
accuracy  is  reported.  Antoon  (Ref  1)  calculated  the 
empirical  variance  of  R(t)  based  on  2000  samples  and  the 
results  in  Table  III  show  empirical  variance  as  a  function 
of  sample  size  and  reliability. 

It  is  useful  to  compare  the  empirical  variance  of 
R(t)  with  the  Cramer-Rao  Lower  Bound  (CRLB).  The  CRLB 
describes  the  lower  limit  for  the  variance  of  an  unbiased 
estimator.  Although  R(t)  is  not  unbiased,  Table  I  does  show 
the  bias  approaching  zero  as  the  sample  size  increases. 

Thus  the  asymptotic  variance  of  R(t)  equals  the  CRLB. 

Rao ' s  results  on  the  asymptotic  distribution  of  a  function 
of  asymptotic  normal  variables  (Ref  10)  can  be  applied  to 
R(t)  . 


CRLB  «a^ 

9 

2 

where a^,  ag  r  ,  and 

6  “  ’ 

covariance  matrix  of 
13:449)  show  that 


+ 


.  aR  aR 
0,k  ae  aK 


2 


are  elements  of  the  asymptotic 
(9,k).  Thoman,  Bain,  and  Antle  (Ref 


=«  1.1099^/(nk^)  -  .608k^/n  a-  C  =  .  257  e/n 

@  iC  0  y  iC 
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VARIANCE  OF  R(t)  (Ref  1) 
(X  10000) 


Therefore  (Ref  14:365) 


CRLB  -  R^  (In  R)  ^  [  1. 109  -  .5141n(-lnR) 

+  .608[ln(-ln  R)  1^] /n  (9) 

One  can  see  that  the  CRLB  approaches  zero  as  the  sample 
size  gets  larger  or  the  reliability  approaches  one.  Snead 
(Ref  11:22)  found  the  CRLB  by  using  the  asymptotic  covari¬ 
ance  matrix  and  numerous  algebraic  operations  involving 

/\  /N 

n,  t,  k,  and  6.  Eq  (9)  represents  the  CRLB  as  a  function 
of  just  R(t)  and  n.  This  is  a  significant  simplification 
over  Snead's  method  without  the  loss  of  any  accuracy.  Eq 
(9)  was  used  to  generate  the  data  in  Table  IV  which  shows 
the  CRLB  as  a  function  of  sample  size  and  reliability.  One 
would  expect  that  the  Table  IV  values  should  be  less  than 
the  corresponding  empirical  variances  in  Table  III.  This 
is  generally  correct,  especially  at  low  reliabilities  and 
small  sample  sizes.  Occasionally,  due  to  Monte  Carlo  vari¬ 
ability,  the  empirical  variance  is  less  than  the  CRLB. 
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III .  Procedure 

The  essential  thrust  of  this  thesis  Is  to  test  a 
univariate  method  of  finding  confidence  intervals  for  the 
reliability  of  complex  systems  over  a  specified  mission 
time.  The  method  will  Incorporate  the  knowledge  gained  In 
Chapter  II  pertaining  to  the  bias  and  empirical  variance 
of  R  (t)  . 


Component  s 

To  fully  test  the  method,  a  variety  of  systems  com¬ 
posed  of  a  variety  of  components  should  be  considered.  The 
true  reliability  of  component  1,  R^,  is  found  by  sub¬ 
stituting  into  the  Weibull  reliability  formula,  Eq  (2). 

In  each  case,  the  mission  time,  t,  is  arbitrarily  set  at 
100  hours.  The  following  components  were  used: 


Component  1 

Failure  Distribution 

Weibull 

Parameter  Values 

k  -  2 

9  -  250 

9  ^ 

■ 

0 

True  Reliability 

Rj^“  exp 

[-(100/250) 

* 

. 85214 

Component  2 

Failure  Distribution 

Weibull 

Parameter  Values 

k  -  3 

9  -  210 

•j  ^ 

a 

0 

True  Reliability 

R^*  exp 

[-(100/210) 

.89765 

Component  3 

Failure  Distribution 

Weibull 

Parameter  Values 

k  -  2 

9  -  300 

9  ^ 

0 

True  Reliability 

R3-  exp 

[-(100/300) 

.89484 
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Component  4 

Failure  Distribution  Welbull 

Parameter  Values  k  ”  3.5  9  »  150  c  *  0 

True  Reliability  ■  exp [ - { 100 / 150) ^ ^ *  .78512 

Component  5 

Failure  Distribution 
Parameter  Values 
True  Reliability 

System  Networks 

The  five  components  were  combined  to  form  various 

kinds  of  systems.  Looking  at  a  variety  of  systems  helps 

to  assure  that  the  results  are  not  unique  to  the  particular 

configuration  of  a  system.  Four  different  systems  were 

used  and  are  shown  in  Fig.  2. 

System  1  consists  of  three  components  in  series. 

The  reliability  of  a  system  with  components  connected  in 

series  is  found  by  taking  the  products  of  the  individual 

component  reliabilities.  Therefore,  the  true  reliability 

of  System  1,  R  ,  is 
®1 

=■  (.85214)  (.89765)  (.89484) 

*  .68448 

System  2  consists  of  one  component  connected  in 
series  with  two  in  parallel.  To  simplify  the  reliability 
equations,  let  be  the  probability  of  failure  for  com¬ 
ponent  i.  Therefore,  »  1-R^.  The  true  reliability  of 


Weibull 

k  -  2.5  9  =  250  c  =«  0 

Rg-  exp[-(100/250)2.5]=  .90376 
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System  1 


System  2 


System  3 


System  4 


SysCem  2,  R  i  Is 

-  .  85214  [  1  -  (.  10235) (.  10516)  ] 

-  .84297 

System  3  consists  of  three  components  connected  in 

parallel.  The  reliability  of  this  system  is  simply  one 

minus  the  probability  that  all  components  fail.  Thus  the 

true  reliability  of  System  3,  R  ,  is 

®3 

^S3  ”  ^  -  ^i‘^2Q3 

=■  1  -  (.14786)  (.10235)  (.10516) 

-  .99841 

All  systems  can  be  seen  as  consisting  of  subsystems 
of  components  which  are  either  series-connected  (such  as 
System  1)  or  parallel-connected  (such  as  System  3).  The 
subsystems  are  in  turn  connected  to  one  another  in  either 
series  or  parallel.  Thus  one  can  combine  subsystem  relia¬ 
bilities  using  either  the  series  formula,  Eq  (10),  or  the 
parallel  formula,  Eq  (12),  until  the  system  reliability 
is  found.  This  was  done  for  System  2  and  will  now  be  done 
for  the  larger  complex  network.  System  4. 
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(13) 


-  R^(l  -  Q2[1  -  Rg(l  -  Q3Q4)]] 

-  .85214[1  -  .  10235  [1  -  .  90376(1  -  ( .  105 16) ( . 2 1488) )  ]  ] 

-  .84197 

The  Monte  Carlo  Method 

Now  the  true  reliabilities  of  all  five  components 
and  all  four  systems  have  been  calculated.  But,  If  the 
reliability  engineer  had  all  the  true  reliabilities,  he 
would  not  need  to  go  through  any  kind  of  Monte  Carlo  simula¬ 
tion.  In  this  study,  the  true  reliabilities  are  considered 
because,  after  the  simulation  is  over,  the  true  reliabili¬ 
ties  will  provide  an  absolute  measure  against  which  the 
univariate  method  can  be  gauged. 

To  find  confidence  intervals  about  the  true  system 
reliabilities,  it  is  desired  to  modify  the  Monte  Carlo 
method  presented  in  Chapter  I  to  incorporate  the  theoretical 
knowledge  gained  in  Chapter  II  about  the  bias  and  empirical 
variance  of  R(t).  The  steps  to  the  Monte  Carlo  method  as 
applied  in  this  thesis  are  as  follows.  A  discussion  of 
each  step  follows  this  list. 

1.  Using  the  true  shape,  scale,  and  location 
parameters,  generate  a  simulated  sample  of  test  data  (com¬ 
ponent  failure  times)  from  the  Weibull  distribution. 

2.  Based  on  the  simulated  test  data,  calculate  the 
maximum  likelihood  estimators  of  the  shape  and  scale 
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parameters.  The  location  parameter  Is  assumed  to  be 
zero  . 

3.  Substitute  the  estimators  of  the  shape  and 
scale  parameters  Into  the  reliability  equation  to  obtain 

a  maximum  likelihood  estimator  of  the  component  reliability. 

4.  Subtract  the  bias  from  the  maximum  likelihood 
estimator  of  the  reliability  to  obtain  an  unbiased  esti¬ 
mator  of  the  component  reliability. 

5.  Form  a  vector  of  component  reliability  esti¬ 
mates  distributed  normally  about  the  unbiased  estimator. 

6.  Repeat  steps  1-5  for  each  component. 

7.  Using  the  reliability  equation  of  a  complex 
system,  combine  the  vectors  of  component  reliabilities  to 
obtain  a  vector  of  system  reliabilities. 

8.  Order  the  vector  of  system  reliabilities  and 
determine  the  99,  95,  90,  80,  70,  60,  and  50  percent  one¬ 
sided  confidence  Intervals.  Note  If  each  of  these  Inter¬ 
vals  contains  the  true  system  reliability. 

9.  Repeat  steps  1-7  until  the  desired  Monte  Carlo 
size  Is  reached. 

10.  To  measure  the  accuracy  of  the  confidence 
limits,  determine  the  percentage  of  the  runs  In  which  each 
of  the  confidence  Intervals  covered  the  true  system  reli¬ 
ability. 

The  computer  program  which  executes  the  above  Monte 
Carlo  method  Is  shown  In  the  appendix.  Some  elaboration  Is 


necessary  to  explain  Individual  steps  of  the  technique. 


In  step  1,  Welbull-dlstributed  sample  failure  times  were 
generated  by  using  the  subroutine  GGWIB  from  the  Inter¬ 
national  Mathematical  and  Statistical  Libraries  (IMSL) . 

/s 

In  step  2,  the  MLEs  of  the  shape  and  scale  parameters  (k 
and  9)  were  determined  by  simultaneous  solution  of  Eq  (5) 
and  Eq  (7)  using  the  nonlinear  simultaneous  equation 
solver  ZSYSTM,  also  from  the  IMSL  library.  In  step  3,  the 
MLE  of  reliability,  R(t),  was  found  by  substitution  of 
k  and  G  Into  Eq  (8).  In  step  4,  the  bias  of  the  MLE  R(t) 
was  determined  by  Interpolation  (using  cubic  splines 
because  of  nonlinearity)  In  Table  I.  Note  that  the  result 
of  step  4  Is  an  estimator  that  Is  usually  less  than  the 
maxlmum-llkellhood  estimator  (because  the  bias  Is  usually 
positive).  This  will  result  in  more  accurate  and  generally 
smaller  component  reliability  estimates  and  therefore 
generally  smaller  estimated  system  reliabilities.  Step  4 
is  an  attempt  to  improve  on  the  optimistic  (high)  system 
reliability  confidence  limits  generated  by  Snead  (Ref  11) 

In  his  univariate  asymptotic  method.  Step  5  Is  a  crucial 
step.  Here  the  assumption  was  made  that  the  unbiased  com¬ 
ponent  reliability  estimates  are  distributed  normally.  The 
mean  Is  assumed  to  be  the  unbiased  R(t)  obtained  from 
step  4.  The  variance  Is  found  by  interpolating  out  of 
Table  III,  again  using  a  cubic  splines  Interpolation 
because  of  nonlinearity  In  the  table.  Since  this  Inter¬ 
polated  variance  should  be  no  less  than  the  Cramer-Rao 
Lower  Bound,  the  CRLB  Is  also  calculated  using  Eq  (9). 
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Then  the  higher  of  the  two  is  used  as  the  variance  to  form 
a  vector  of  2000  random  and  normally-distributed  component 
reliability  estimates.  (The  IMSL  subroutine  GGNML  was  used 
to  generate  random  normal  deviates.)  In  step  7,  the 
vectors  of  component  reliability  estimates  are  combined 
using  the  systems  reliability  equations  [Eqs  (10),  (11), 
(12),  and  13)].  The  result  is  a  vector  of  2000  estimated 
system  reliabilities  for  each  system.  These  vectors  are 
each  ordered  in  ascending  sequence  and  then  the  1,  5,  10, 


20, 

30, 

40. 

and 

50 

percent iles 

are  chosen  to 

determine  the 

99. 

95, 

90, 

80  , 

70 

,  60,  and  50 

percent  lower 

confidence 

limits.  Each  of  these  lower  confidence  limits  is  compared 
to  the  true  system  reliability.  If  the  lower  confidence 
limit  is  less  than  or  equal  to  the  true  system  reliability, 
then  the  associated  confidence  interval  contains  the  true 
system  reliability. 

The  entire  experiment  is  repeated  again  and  again 
(step  9)  until  the  desired  Monte  Carlo  size  is  reached. 

In  this  study,  600  runs  of  the  simulation  were  made.  In 
each  run,  and  for  each  system  and  each  confidence  interval, 
it  was  noted  if  the  confidence  interval  contained  the  true 
system  reliability.  Ideally,  the  X  percent  confidence 
interval  should  contain  the  true  system  reliability  X  per¬ 
cent  of  the  time.  For  example,  it  is  desired  that  the  95 
percent  confidence  interval  contain  the  true  system  reli¬ 
ability  during  95  percent  of  the  600  runs.  This  provides 
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an  absolute  test  (step  10)  of  the  accuracy  of  this  unl~ 
variate  method. 


Additional  Testing  and 
Senslt Iv Ity  Analysis 

In  an  attempt  to  validate  the  results  of  the  Monte 
Carlo  technique,  some  additional  testing  was  conducted.  One 
test  was  to  check  the  assumption  of  normality  made  In 
step  5.  It  Is  already  known  that  estimates  of  component 
reliability  are  only  asymptotically  normal.  A  Kolmogorov- 
Smlrnoff  test  was  run  to  check  the  distribution  of  com¬ 
ponent  reliability  estimates. 

Even  though  Table  II,  the  variances  of  R(t)  given 
by  Thoman,  Bain,  and  Antle  (Ref  14:366),  contains  only 
three  significant  digits  at  most.  It  was  decided  to  check 
If  the  results  were  sensitive  to  the  use  of  that  table. 

Also,  since  a  faulty  random  number  generator  can  signifi¬ 
cantly  affect  the  results  of  any  slmultation,  a  random 
normal  generator  other  than  GGNML  was  tried.  GGNML  uses 
the  Inverse  normal  function  of  uniform  deviates.  The  Box- 
Muller  (Ref  3:610-11)  was  also  utilized.  This  technique 
uses  the  natural  logarithm  and  trigonometric  functions  to 
generate  normal  deviates  from  random  uniform  variables. 

Since  GGNML  and  Box-Muller  use  different  uniform  deviate 
generators,  the  generation  of  uniform  deviates  Is  also 
being  checked  by  two  methods.  Another  desirable 
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characteristic  of  a  simulation 
to  the  choice  of  the  initial  r 
simulations  were  run  with  dlff 
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IV.  Results 


Base  Case 

Tables  V  through  VIII  contain  the  base  case  results 
of  the  Monte  Carlo  simulation  for  Systems  1  through  4. 

For  each  system,  the  simulation  was  run  for  component  test 
data  sample  sizes  of  10,  15,  20,  50,  and  100.  Each  table 
entry  is  the  percentage  of  the  Monte  Carlo  runs  in  which 
the  simulated  confidence  interval  covered  the  true  system 
reliability.  For  example,  in  the  first  row  of  Table  V  and 
with  a  sample  size  of  10  failure  times  for  each  component, 
the  99  percent  confidence  interval  contained  the  true  reli¬ 
ability  of  System  1  only  .9383  of  the  runs.  As  the  com¬ 
ponent  sample  size  increases  in  that  first  row,  one  can 
see  that  the  confidence  interval  coverage  improves  until, 
at  a  sample  size  of  100,  the  coverage  is  .9817,  almost  99 
percent.  This  general  improvement  is  also  the  case  when 
looking  at  increasing  sample  sizes  for  the  other  confidence 
intervals . 

Tables  VI  through  VIII  reveal  that  Systems  2,  3, 
and  4  experience  the  same  improvements  as  component  sample 
sizes  get  larger.  For  each  system,  the  confidence  interval 
coverage  is  low  for  small  sample  sizes.  The  reason  for  this 
is  high  or  optimistic  system  reliability  estimates.  As  a 
result,  the  lower  confidence  limits  are  also  too  high  and 
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TABLE  V 


SYSTEM  1  CONFIDENCE  INTERVAL  COVERAGE  OF  THE 
TRUE  SYSTEM  RELIABILITY 


Confidence  Interval 
(Percent ) 

10 

15 

Sample  Size 
20 

50 

100 

99 

.  9383 

.  9450 

.  9617 

.  9733 

.  9817 

95 

.  8867 

.8967 

.  9050 

.  9267 

.  9333 

90 

.  8267 

.  8633 

.  8617 

.  8783 

.  8950 

80 

.  7417 

.  7683 

.  7800 

.  7983 

.  8233 

70 

.6650 

.6683 

.6817 

.7183 

.  7200 

60 

.  5750 

.5783 

.  5700 

.6167 

.6367 

50 

.4917 

.5083 

.4833 

.5067 

.  5333 

TABLE  VI 

SYSTEM  2  CONFIDENCE  INTERVAL  COVERAGE  OF  THE 
TRUE  SYSTEM  RELIABILITY 


Confidence  Interval 
(Percent) 

10 

15 

Sample  Size 
20 

50 

100 

99 

.8717 

.  9100 

.  9383 

.  9450 

.  9667 

95 

.  7933 

.8467 

.  8600 

.  9083 

.  9200 

90 

.  7467 

.7917 

.8083 

.  8617 

.  8817 

80 

.6750 

.7017 

.7200 

.  7767 

.7917 

70 

.6050 

.6033 

.6467 

.6750 

.7117 

60 

.5083 

.5067 

.  5433 

.5567 

.  5950 

50 

.4283 

.4283 

.  4700 

.  4583 

.5067 
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TABLE  VII 

SYSTEM  3  CONFIDENCE  INTERVAL  COVERAGE  OF  THE 
TRUE  SYSTEM  RELIABILITY 


Confidence  Interval 
(Percent) 

10 

15 

Sample  Size 
20 

50 

100 

99 

.  7967 

.  8717 

.  9017 

.  9467 

.  9717 

95 

.6517 

.  7867 

.  7950 

.8700 

.  9017 

90 

.5917 

.  7033 

.  7050 

.8100 

.  8467 

80 

.4633 

.5283 

.  5733 

.6850 

.  7467 

70 

.  3600 

.4183 

.  4617 

.  5783 

.6133 

60 

.2667 

.  3200 

.  3650 

.4683 

.5200 

50 

.  1850 

.2267 

.  2767 

.  3683 

.  4550 

SYSTEM  4 

TABLE  VIII 

CONFIDENCE  INTERVAL  COVERAGE 
TRUE  SYSTEM  RELIABILITY 

OF  THE 

Confidence  Interval 

(Percent)  10 

Sample  Size 

15  20  50 

100 

99 

.  8750 

.  9100 

.  9417 

.  9483 

.  9750 

95 

.  7983 

.8433 

.  8567 

.9067 

.  9200 

90 

.7400 

.7933 

.  8033 

.  8583 

.  8817 

80 

.6783 

.7067 

.7183 

.  7733 

.  7767 

70 

.5967 

.6100 

.  6450 

.  6633 

.  7067 

60 

.5100 

.5083 

.5533 

.5583 

.5983 

50 

.4350 

.4233 

.  4733 

.4550 

.4950 
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the  confidence  Interval  Is  not  wide  enough  to  cover  the 
true  system  reliability. 

An  unexpected  result  from  the  simulation  was  the 
effect  of  true  system  reliability  on  the  confidence  inter¬ 
val  coverage.  It  appears  that  a  lower  system  reliability 
corresponds  to  a  more  accurate  confidence  interval. 

System  1,  with  a  reliability  of  .68448,  has  consistently 
more  accurate  confidence  interval  coverage  than  any  other 
system.  Systems  2  and  4  both  have  reliabilities  of  about 
.84  and  corresponding  entries  in  Table  VI  and  VIII  are 
very  close.  System  3  has  the  highest  reliability,  .99841, 
and  also  the  least  accurate  confidence  intervals.  Even  at 
a  sample  size  of  100,  Table  VII  shows  the  60  and  70  percent 
confidence  interval  coverage  being  eight  to  nine  percent 
too  low.  Other  systems  have  confidence  intervals  for  a 
sample  size  of  100  that  are  always  less  than  four  percent 
away  from  the  ideal  accuracy. 

Sensitivity  Analysis 

Sensitivity  analysis  was  conducted  to  test  some  of 
the  model  Inputs.  To  build  confidence  in  the  results  of  a 
simulation,  it  is  desired  that  the  results  not  be  sensi¬ 
tive  to  the  choice  of  random  number  streams.  Therefore, 
a  different  initial  random  number  seed  was  chosen  and  the 
results  for  all  four  systems  are  shown  in  Table  IX.  The 
95,  90,  and  80  percent  confidence  intervals  were  calculated 
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TABLE  IX 

CONFIDENCE  INTERVAL  COVERAGE  WITH  A 
DIFFERENT  RANDOM  NUMBER  STREAM 


Confidence  Interval 
(Percent) 

10 

15 

Sample  Size 
20 

50 

100 

System  1 

95 

.8917 

.  8883 

.  9217 

.  9233 

.  9467 

90 

.8350 

.8317 

.  8917 

.8717 

.  8850 

80 

.7433 

.7417 

.7833 

.  7783 

.7850 

System  2 

95 

.8300 

.  8533 

.  8500 

.  8900 

.9300 

90 

.  7800 

,  8050 

.  8017 

.  8417 

.  8650 

80 

.6883 

.6917 

.  7067 

.7500 

.  7750 

System  3 

95 

.  6750 

.7433 

.  8133 

.  8550 

.  9050 

90 

.6050 

.6683 

.  7250 

.  7867 

.8300 

80 

.4650 

.5200 

.5833 

.6667 

.7150 

System  4 

95 

.  8233 

.8583 

.  8517 

.8867 

.  9233 

90 

.7817 

.8067 

.  7983 

.8383 

.8683 

80 

.6867 

.7017 

.  7050 

.7533 

.7717 

1 


« 
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and  the  table  reveals  only  slight  variations  from  the  cor¬ 
responding  entries  in  Tables  V  through  VIII. 

A  faulty  random  number  generator,  especially  one 
that  does  not  generate  enough  extreme  values  of  the  dis¬ 
tribution,  can  significantly  affect  the  results  of  this 
simulation.  Thus  Che  IMSL  random  normal  generator  GGNML 
was  replaced  with  the  Box-Muller  normal  generator.  The 
Box-Muller  method  is  exact.  That  is,  if  the  uniform 
generator  is  truly  generating  random  uniform  deviates, 
then  the  Box-Muller  equations  (Ref  3:610)  will  output  truly 
random  normal  deviates.  This  will  guarantee  a  proper  dis¬ 
tribution  in  the  tails  of  the  normal  density  function. 
Simulation  runs  at  sample  sizes  of  10  and  15  were  made  and 
confidence  interval  coverage  is  shown  in  Table  X.  Again, 
only  small  variations  (due  to  Monte  Carlo  variability) 
exist  between  this  table  and  Tables  V  through  VIII. 

Thoman,  Bain,  and  Antle  calculated  the  empirical 
variances  displayed  in  Table  II.  Although  they  reported 
only  three  significant  digits,  they  did  use  a  Monte  Carlo 
size  of  10,000  as  opposed  to  Table  III,  in  which  more  sig¬ 
nificant  digits  are  shown,  but  the  Monte  Carlo  size  was 
2000.  Thoman,  Bain,  and  Antle’s  table  was  used  to  see  if 
significant  variations  would  occur.  The  95,  90,  and  80 
percent  confidence  Interval  coverage  is  shown  in  Table  XI 
for  sample  sizes  of  10,  15,  20,  and  50.  Differences 
between  these  results  and  the  base  case  results  are  minor. 
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TABLE  X 


CONFIDENCE  INTERVAL  COVERAGE  WITH  THE  BOX-MULLER 
RANDOM  NUMBER  GENERATOR 


Confidence  Confidence 


Interval 

(Percent) 

S  ample 
10 

Size 

15 

Interval 

(Percent) 

S  amp le 
10 

Size 

15 

System  1 

System  3 

99 

.  9433 

.  9467 

99 

.  7950 

.8667 

95 

.  8783 

.8933 

95 

.6733 

.7367 

90 

.  8217 

.  8300 

90 

.5900 

.6650 

80 

.  7367 

.  7433 

80 

.4417 

.5117 

70 

.6517 

.6533 

70 

.  3667 

.4250 

60 

.5533 

.5733 

60 

.2783 

.3033 

50 

.  4850 

.4750 

50 

.1883 

.2183 

System  2 

System  4 

99 

.  8650 

.  9050 

99 

.8583 

.  9083 

95 

.  8150 

.8333 

95 

.  8167 

.  8383 

90 

.7533 

.  7750 

90 

.  7633 

.7767 

80 

.6867 

.  7000 

80 

.6900 

.6967 

70 

.6017 

.5900 

70 

.  5950 

.5967 

60 

.5017 

.5267 

60 

.  5000 

.5183 

50 

.4100 

.4283 

50 

.4033 

.  4367 

TABLE  XI 

CONFIDENCE  INTERVAL  COVERAGE  WITH 
VARIANCES  OF  THOMAN ,  BAIN,  AND 

EMPIRICAL 

ANTLE 

- 

Confidence  Interval 

Sample 

Size 

(Percent) 

10 

15 

20 

50 

System  1 

i 

I 

f 

95 

.8917 

.  8817 

.9217 

.  9233 

1 

90 

.8350 

.  8300 

.8900 

.  8750 

80 

.7417 

.7400 

.  7833 

.7800 

System  2 

95 

.  8300 

.  8533 

.  8500 

.8900 

90 

.  7800 

.  8050 

.  8017 

.  8417 

80 

.6883 

.6867 

.7067 

.7500 

System  3 

95 

.6750 

.  7417 

.  8133 

.  8583 

90 

.6050 

.6650 

.  7250 

.7867 

80 

.4650 

.5200 

.5817 

.6667 

, 

System  4 

95 

.  8233 

.  8583 

.8517 

.  8867 

90 

.7817 

.8083 

.7983 

.8383 

' 

1 

80 

.6867 

.6950 

.  7050 

.  7533 

i 

i 

\ 

> 

t 

1 

1 

•  : 

■ 

: 

; 
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The  sensitivity  analysis  has  revealed  only  small 
variations  (a  maximum  of  less  than  five  percent)  between 
the  base  case  and  choice  of  random  number  stream,  type  of 
normal  and  uniform  random  number  generator,  and  choice  of 
empirical  variance  table.  These  slight  differences  can 
occur  simply  as  a  result  of  the  variability  of  Monte  Carlo 
simulation  (as  In  any  random  process) .  A  more  fundamental 
question  Is  the  assumption  of  normality  made  In  step  5  of 
the  Monte  Carlo  method. 

Ko Imo go rov- Smirnoff  Test 
of  Normality 

Since  the  Monte  Carlo  size  was  600,  there  were 
also  600  samples  of  failure  times  drawn  for  each  of  the 
five  components.  In  turn,  the  sample  sizes  were  either  10, 
15,  20,  50,  or  100.  As  a  result,  for  each  component  and 
for  each  of  the  sample  sizes,  there  were  600  estimates  of 
the  component  reliability.  A  Ko Imogorov-Smirnof f  test  was 
run  for  the  distribution  of  these  reliability  estimates. 

The  null  hypothesis  was  that  the  estimates  are  normally  dis¬ 
tributed  with  a  mean  equal  to  the  true  component  reliabil¬ 
ity  and  a  variance  equal  to  the  interpolated  value  from 
Table  III  or  the  Cramer-Rao  Lower  Bound  (whichever  Is 
larger).  The  Ko Imogorov-Smirnof f  test  was  made  with  a  .05 
probability  of  Type  I  error.  That  is,  if  the  distribution 
Is  actually  from  the  normal  density  function,  the  test 
still  has  a  .05  probability  of  rejecting  the  null  hypothesis. 
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Table  XII  contains  the  results  of  the  test  as  well  as  the 


bias  of  the  estimated  reliability  and  the  difference  between 
the  sample  standard  deviation  and  the  assumed  population 
standard  deviation.  There  Is  a  tendency  to  reject  normality 
for  sample  sizes  of  20  or  less  and  to  accept  normality  for 
sample  sizes  of  50  or  greater.  The  table  also  suggests 
that  the  reason  for  rejecting  normality  Is  skewness  In  the 
distribution  of  R(t)  for  small  sample  sizes.  For  example. 
Component  3  was  rejected  at  a  sample  size  of  15  (with  a 
bias  of  .00132  and  a  difference  of  standard  deviations  of 
-.00011),  while  Component  2  passed  the  normality  test  at 
a  sample  size  of  100  (with  a  bias  of  -.00165  and  a  differ¬ 
ence  In  standard  deviations  of  .00139).  If  the  distribu¬ 
tion  of  R(t)  Is  skewed  for  small  sample  sizes,  rejection 
of  normality  would  still  occur  even  though  the  bias  and 
difference  In  standard  deviations  Is  so  small.  There  Is  a 
reasonable  explanation  for  this  skewness.  Since  the  reli¬ 
ability  estimates  are  bounded  by  zero  and  one,  as  the  true 
component  reliability  approaches  one,  the  distribution  of 
the  reliability  estimates  Is  significantly  affected  by  the 
upper  bound  of  one,  but  not  by  the  lower  bound  of  zero. 
Therefore,  the  distribution  would  not  be  symmetric,  but  be 
skewed  to  the  left. 
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TABLE  XII 


KOLMOGOROV-SMIRNOFF  TEST  OF  NORMALITY  FOR  DISTRIBUTION 
OF  COMPONENT  RELIABILITY  ESTIMATES 


10 

15 

Sample  Size 

20  50 

100 

Component  1 

Reject /Pass 

Reject 

Rej  ect 

Reject 

Pass 

Pass 

E(R)  -  R 

.00276 

.00179 

.00072 

.00031 

-.00003 

(s  -  a) 

.00439 

.00133 

-.00034 

.00149 

-.00016 

Component  2 

ReJ  ect /Pass 

Reject 

Reject 

Reject 

Pass 

Pass 

E(R)  -  T 

.00199 

-.00227 

-.00131 

-.00073 

-.00165 

(s  -  a) 

-.00120 

.00303 

.00221 

.00072 

.00139 

Component  3 

Rej  ect /Pass 

Reject 

Reject 

Rej  ect 

Pass 

Pass 

E(R)  -  R 

.00296 

.00132 

.00077 

-.00117 

-.00086 

(s  -  a) 

-.00118 

-.00011 

.00209 

-.00137 

.00018 

Component  4 

Rej  ect /Pass 

Reject 

Reject 

Pass 

Pass 

Pass 

E(R)  -  R 

.00339 

.00834 

.00278 

-.00284 

.00087 

(s  -  a) 

.00068 

-.00266 

-.00352 

-.00024 

-.00045 

Component  5 

Reject/Pass 

Reject 

Reject 

Reject 

Pass 

Rej  ect 

E(R)  -  R 

.00121 

-.00137 

.00239 

-.00031 

.00120 

(s  -  CT) 

.00044 

.00237 

.00135 

.00011 

-.00019 
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V. 


Cone lus ions  and  Recommendations 


Conclus ions 

The  univariate  method  developed  in  this  thesis 
works  best  when  component  reliabilities  and  system  relia¬ 
bilities  are  not  high  (say,  less  than  .9).  If  component 
reliabilities  are  higher  than  .9,  then  the  skewness  of  the 

A 

distribution  of  R(t)  has  a  definite  impact  on  the  results. 
If  the  system  reliability  is  greater  than  .9,  then  it  is 
very  sensitive  to  even  slight  changes  in  component  reliabil 
ity  estimates.  Therefore,  even  a  slight  skewness  will  sig¬ 
nificantly  affect  the  system  reliability  estimates.  If 
high  component  or  system  reliabilities  exist,  then  they  can 
be  compensated  for  by  large  sample  sizes.  A  sample  size  of 
50  or  greater  is  adequate  unless  the  component  or  system 
reliabilities  exist,  then  they  can  be  compensated  for  by 
large  sample  sizes.  A  sample  size  of  50  or  greater  is 
adequate  unless  the  component  or  system  reliabilities  are 
extremely  close  to  one. 

Recommendations 

The  most  immediate  suggestion  for  further  work  is 
that  investigation  be  conducted  on  the  distribution  of  com- 

A 

ponent  reliability  estimates,  R(t),  based  on  sample  sizes 
of  less  than  50.  The  distribution  is  definitely  not  normal 
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Therefore,  it  would  be  helpful  to  try  fitting  other  dis- 

A 

tributions  to  the  empirical  distribution  of  R(t).  Since 
(t)  is  bounded  between  zero  and  one,  a  logical  choice  is 
the  Beta  distribution.  But  most  importantly,  the  Beta  can 
also  fit  the  skewness  of  the  empirical  distribution. 

Because  many  critical  components  and  systems  are  built  to 
demonstrate  extremely  high  reliabilities,  particular  empha¬ 
sis  should  be  placed  on  component  reliabilities  of  greater 
than  . 9 . 

Once  a  close  fit  is  found  for  the  distribution  of 
component  reliability  estimates,  the  next  step  is  Monte 
Carlo  simulation  of  the  reliability  of  entire  systems.  The 
computer  program  shown  in  the  appendix  has  the  capability 
to  execute  the  systems  simulation.  Only  slight  modifica¬ 
tions  are  required  to  change  the  sampling  distribution  (say 
the  Beta  Instead  of  the  normal) ,  component  test  data  sample 
size,  true  component  reliabilities,  or  the  configuration 
of  systems. 


44 


1. 


Bibliography 


Antoon,  David  F.  Conf idence  Intervals  and  Tests  of 
Hypothesis  of  £  System  whose  Underlying  Distribution 
is  ^  Two  Parameter  Weibull.  Unpublished  thesis.  Wright- 
Patterson  Air  Force  Base,  Ohio:  Air  Force  Institute  of 
Technology,  December  1979. 

2.  Bernhoff,  Albert  0.  Conf idence  Limits  for  System 
Reliability  Based  on  Component  Test  Data .  Unpublished 
thesis.  Wr ight-Patter son  Air  Force  Base,  Ohio:  Air 
Force  Institute  of  Technology,  August  1963. 

3.  Box,  G.  E.  P.  and  Mervin  E.  Muller.  "A  Note  on  the 
Generation  of  Random  Normal  Deviates,"  Annals  of 
Mathematical  Statistics  2_9:610,11  (1958). 

4.  Lannon,  Robert  G.  A  Mont e  Carlo  Technig ue  for  Approxi¬ 
mating  System  Reliability  Confidence  Limits  Using  the 
Weibull  Distribution .  Unpublished  thesis.  Wright- 
Patterson  Air  Force  Base,  Ohio:  Air  Force  Institute  of 
Technology,  1972. 

5.  Levy,  Louis  L,  and  Albert  H.  Moore.  "A  Monte  Carlo 
Technique  for  Obtaining  System  Reliability  Confidence 
Limits  from  Component  Test  Data,"  IEEE  Transactions  on 
Reliability ,  R-16,  ^:69-72  (September  1967). 

6.  Mendenhall,  William  and  Richard  L.  Scheaffer.  Mathe¬ 
matical  Statistics ,  with  Applications .  North  Scituate, 
Massachusetts:  Duxbury  Press,  1973. 

7.  Moore,  Albert  H.  "Extension  of  Monte  Carlo  Technique 
for  Obtaining  System  Reliability  Confidence  Limits  from 
Component  Test  Data,"  Proceedings  of  National  Aerospace 
Electronics  Conference ,  459-463  (May  1965)  . 

8.  Orkand,  Donald  S.  A  Monte  Carlo  Method  for  Determining 
Lower  Confidence  Limits  for  System  Reliability  on  the 
Basis  of  Sample  Component  Data .  Report  No.  ORDBB-VC-4. 
Dover,  New  Jersey;  Concepts  and  Applications  Laboratory, 
Picatinny  Arsenal,  June  1960.  (AD  627  799). 

9.  Park,  Won  J.  Basic  Concepts  of  Statist ics  and  their 
Applications  in  Composite  Materials .  Technical  Report. 
Wright-Patterson  Air  Force  Base,  Ohio:  Air  Force 
Materials  Laboratory,  1979. 


45 


10 


I 


Methods  In 
Sons ,  1952 . 

11.  Snead,  Robert  C.  A  Univariate  Asymptotic  Normal  Monte 
Carlo  Method  (UANMCM)  for  Estimating  System  Reliability 
Confidence  Limits .  Unpublished  thesis.  Wrlght- 
Patterson  Air  Force  Base,  Ohio:  Air  Force  Institute  of 
Technology,  December  1978. 

12.  T echnlcal  Tralnlng--Rel lab 11 Ity /Maintainability . 

Chanute  Air  Force  Base,  Illinois:  Chanute  Technical 
Training  Center,  March  1978. 


Rao ,  C.  Radhakr ishna .  Advanced  Statistical 
Biometric  Research .  New  York:  John  Wiley  & 


13.  Thoman,  Darrel  R.,  Lee  J.  Bain  and  Charles  E.  Antle. 
"Inferences  on  the  Parameters  of  the  Weibull  Distribu¬ 
tion,"  Technometrics ,  lA:445-60  (August  1969). 

14.  Thoman,  Darrel  R.,  Lee  J.  Bain  and  Charles  E.  Antle. 
"Maximum  Likelihood  Estimation,  Exact  Confidence  Inter¬ 
vals  for  Reliability,  and  Tolerance  Limits  in  Weibull 
Distribution,"  Tec hnome tries,  1^:363-71  (May  1970). 

15.  Zelen,  M.  and  M.  C.  Dannemiller.  "The  Robustness  of 
Life  Testing  Procedures  Derived  from  the  Exponential 
Distribution,"  Technometrics .  ^:29-49  (February  1961). 


46 


Appendix 

Computer  Program  Listing 


PROGRAM  WEIB  (TAPEl,  INPUT,  OUTPUT) 

C 

c  *  * 

C  *  FOR  EACH  COMPONENT  OF  A  COMPLEX  SYSTEM,  THIS  PROGRAM  * 

C  *  GENERATES  A  SAMPLE  OF  SIZE  NSAM  FROM  THE  WEIBULL  DISTRIBUTION  * 

C  *  USING  THE  TRUE  PARAMETERS  TK  (SHAPE),  TTHETA  (SCALE),  AND  TC  * 

C  *  (LOCATION).  FROM  THIS  SAMPLE,  THE  MAXIMUM  LIKELIHOOD  * 

C  *  ESTIMATORS  (MLES)  FOR  K  AND  THETA  ARE  DERIVED  USING  HARTER  &  * 

C  *  MOORES  ITERATIVE  SCHEME  AND  THE  SIMULTANEOUS  EQUATION  SOLVER  * 

C  *  ZSYSTM.  THE  MLES  ARE  COMBINED  TO  YIELD  RCHAT,  THE  MLE  FOR  THE  * 

C  *  COMPONENT  RELIABILITY.  GIVEN  THE  RELIABILITY  AND  THE  SAMPLE  * 

C  *  SIZE  NSAM,  RCHAT  IS  ASYMPTOTICALLY  NORMALLY  DISTRIBUTED  WITH  A  * 

C  *  SPECIFIED  VARIANCE.  THEREFORE  WE  CAN  SAMPLE  FROM  THE  NORMAL  * 


C  *  DISTRIBUTION  TO  OBTAIN  A  VECTOR  OF  SAMPLE  COMPONENT  * 
C  *  RELIABILITIES.  THIS  PROCESS  IS  REPEATED  FOR  EACH  COMPONENT  * 
C  *  AND  THEN  THE  RELIABILITIES  ARE  COMBINED  FOR  4  DIFFERENT  TYPES  * 
C  *  OF  SYSTEMS  TO  YIELD  4  VECTORS  OF  SAMPLE  SYSTEM  RELIABILITES.  * 


C  *  THESE  VECTORS  ARE  ORDERED  AND  THEN  THE  95,  90,  AND  80  PERCENT  * 

C  *  LOWER  CONFIDENCE  LIMITS  ARE  PICKED.  THIS  ESTABLISHES  THE  95,  * 

C  *  90,  AND  80  PERCENT  CONFIDENCE  INTERVALS  FOR  EACH  SYSTEM  AND  IT  * 

C  *  IS  NOTED  WHETHER  EACH  OF  THESE  INTERVALS  CONTAINS  THE  TRUE  * 

C  *  SYSTEM  RELIABILITY.  * 

C  *  THE  ABOVE  PROCESS  IS  REPEATED  FOR  NOLMC  MONTE  CARLO  RUNS,  WITH  * 

C  *  COUNTERS  FOR  EACH  SYSTEM  TO  TRACK  THE  NUMBER  OF  TIMES  THAT  THE  * 

C  *  THE  ABOVE  PROCESS  IS  REPEATED  FOR  NOLMC  MONTE  CARLO  RUNS,  WITH  * 

C  *  COUNTERS  FOR  EACH  SYSTEM  TO  TRACK  THE  NUMBER  OF  TIMES  THAT  THE  * 

C  *  CONFIDENCE  INTERVALS  CONTAIN  THE  TRUE  SYSTEM  RELIABILITY.  * 

C  *  * 

c 

DIMENSION  ARCHAT  (2000,5),  ATRS  (4),  BI  (12,15),  C  (4,4),  DEV 

1  (2000),  PARAM  (3,5),  QNTKS  (40),  R  (205,5),  RC  (2000,5),  RCHAT 

2  (5),  RCHSIG  (5),  RLBS  (12),  RS  (2000),  SIG  (12,15),  SIGSQ 

3  (12,15),  SMSZ  (15),  TEMP  (200),  TRC  (5),  TRS  (1),  VRH  (10,11), 

4  WK  (200),  WKAREA  (10),  X  (2) 

EXTERNAL  F 

DOUBLE  PRECISIONDSEED 

REAL  KHAT 

DATA  RLBS  /  .5,  .55,  .6,  .65,  .7,  .75,  .8,  .85,  .9,  .925,  .95, 

1  .98  / 

C 

DATA  SMSZ  /  8.,  9.,  10.,  11.,  12.,  13.,  14.,  15.,  20.,  25., 

1  30.,  40.,  50.,  75.,  100.  / 

C 


47 


mi 


o  o 


DATA 

VRH 

/  266., 

242. 

,  194. 

,  163 

• » 

130., 

95 

•  9 

59., 

41., 

25., 

1 

6., 

200. , 

187. 

.  153., 

130., 

103., 

76 

•  » 

47. 

9 

33. 

19 

. ,  5. , 

167 

2 

154. 

,  126 

. ,  107., 

86., 

63., 

39., 

27. 

f 

16., 

4 

•  t 

124. 

,  118. 

.  99 

3 

86. , 

70. , 

51., 

32. 

.  22 

. ,  13. 

.  3., 

90 

•  9 

86. 

9 

72. 

.,  62 

.,  51. 

.  37 

4 

23., 

16. , 

9., 

2., 

72., 

68., 

58., 

50. 

» 

41., 

30. , 

,  19. 

,  13., 

7., 

5 

2.. 

59., 

57., 

48., 

42. 

.  34., 

25., 

16 

•  9 

11. 

9 

6., 

.  1., 

43., 

42., 

6 

36., 

31., 

26., 

19. 

,  12 

. ,  8. , 

5., 

1., 

34. , 

33 

•  9 

29., 

25., 

20. , 

7 

15., 

9., 

7.,  4 

•.  1 

.,  23.,  22 

.,  19 

• , 

17 

.,  14. 

,  10., 

6.,  4. 

,  3. 

8 

1., 

17., 

16., 

14., 

12. 

,  10., 

8., 

5., 

3 

2 

•  9 

1. 

.  / 

DATA 

QNTKS  / 

.975 

,  .842,  .708,  . 

624 

9 

.563 

9 

.5] 

L9,  . 

483,  . 

454, 

1 

.430, 

.409, 

.391, 

.375, 

.361, 

.349, 

.338, 

.327, 

.318, 

.309, 

2 

.301, 

.294, 

.287, 

.281, 

.275, 

.269, 

.264, 

.259, 

.254, 

.250, 

3 

.246, 

.242, 

.238, 

.234, 

.231, 

.227, 

.224, 

.221, 

•> 

00 

CM 

• 

.215, 

4 

.213, 

.210 

/ 

o  o  o  o 


NS4C60 

NS4C50 


-  0 
-  0 


C 

PRINT  409 

C  READ  THE  SAMPLE  SIZE  AND  NUMBER  OF  MONTE  CARLO  RUNS 
READ  *,  NSAM,  NOLMC 
PRINT  419,  NSAM,  NOLMC 
RNSAM  -  NSAM 

C  READ  THE  TRUE  COMPONENT  PARAMETERS 

READ  *,  ((PARAMd,  J) ,  I-l,  3),  J-1,  5) 

PRINT  429,  ((J,  (PARAMd,  J) ,  I-l,  3),  CMPREL(TIME,  PARAMd,  J), 
1  PARAM(2,  J),  PARAM(3,  J))),  J-1,  5) 

READ  (1,  *)  ((SIGd,  J),  J-1,  15),  1=1,  12) 

READ  (1,  *)  ((Bid,  J),  J-1,  15),  I-l,  12) 

DO  10  I  -  1,  12 
DO  10  J  -  1,  15 

SIGSQd,  J)  -  SIGd,  J)  *  SIGd,  J)  *  1.E4 
10  CONTINUE 

PRINT  439,  ((SIGd,  J),  J-1,  15),  I-l,  12) 

PRINT  439,  ((Bid,  J),  J-1,  15),  1=1,  12) 

C 

Q  AAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAA AAAAAAAAAAAAAAAAA AAA AAAA AAAA AA AAAA 

C  THIS  OUTSIDE  LOOP  FROM  HERE  TO  STATEMENT  300  COMPLETES 
C  NOLMC  MONTE  CARLO  RUNS  OF  THE  SIMULATION. 

C 

DO  300  NCOUNT  -  1,  NOLMC 


FOR  EACH  OF  5  COMPONENTS,  THIS  LOOP  GENERATES  THE  MLE  OF  RCHAT 
C  FROM  THE  MLES  OF  K  AND  THETA.  IT  DETERMINES  THE  VARIANCE  OF 
C  RCHAT  AND  THEN  SAMPLES  FROM  THE  NORMAL  DISTRIBUTION  FOR  NDEV 
C  SAMPLE  RELIABILITIES  OF  EACH  COMPONENT. 

C 


DO  290  J 

-  1.  5 

TK 

=  PARAMd,  J) 

TTHETA 

-  PARAM(2,  J) 

TC 

-  PARAM(3,  J) 

C  PASS  INFORMATION  ON  COMPONENT  NUMBER  AND 
C  SAMPLE  SIZE  THRU  ARRAY  R 
R(205,  1)  -  J 
R(204,  J)  -  NSAM 

C  DETERMINE  THE  TRUE  COMPONENT  RELIABILITY 

TRC(J)  -  CMPREL(TIME,  TK,  TTHETA,  TC) 
C 

C  GATHER  NSAM  FAILURE  TIMES  FROM  THE 
C  WEIBULL  DISTRIBUTION  WITH  TRUE  PARAMETERS 
C  TK, TTHETA,  AND  TC 
C 

CALL  GGWIB  (DSEED,  TK,  NSAM,  TEMP) 

DO  20  I  -  1,  NSAM 
R(I,  J)  -  TTHETA  *  TEMP(I)  +  TC 
20  CONTINUE 


c 

c 

c 

c 

c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


DETERMINE  THE  LIKELIHOOD  OF  DRAWING  THIS  SAMPLE 
USING  THE  TRUE  K,  THETA,  AND  C 

DETERMINE  THE  MLE  OF  THETA  AND  K 
BY  HARTER  &  MOORES  ITERATIVE  SCHEME 

(ZSYSTM  SOLVES  THE  TTro  NONLINEAR  SIMULTANEOUS  EQUATIONS) 


30 


40 


EPS 

NSIG 

N 

NSTART 

START 

ITMAX 

SUMXIK 

X(l) 

DO  40  I 
SUMXIK 
X(2) 

CALL  ZSYSTM 

KHAT 

THAT 


l.E  -  10 
8 
2 
0 

2.55 

50 

0. 

START 

-  1.  NSAM 

-  SUMXIK  +  R(I,  J)  **  X(l) 

(SUMXIK  /  NSAM)  **  (1.  /  X(l)) 

(F,  EPS,  NSIG,  N,  X,  ITMAX,  WKAREA, 
X(l) 

X(2) 


R, 


lER) 


IF  (ITMAX  .EQ.  50  .OR.  KHAT  .LT.  0.  .OR.  THAT  .LT.  0.  .OR. 
1  lER.NE.O)  GO  TO  200 

AND  PRINT  ERROR  MESSAGE  THEN  RESTART  ITERATION 
WITH  A  DIFFERENT  KHAT 


DETERMINE  RHAI,  THE  MLE  OF  THE  COMPONENT  RELIABILITY,  USING 
THE  TRUE  C  AND  THE  MLES  OF  K  AND  THETA 

RCHAI(J)  -  CMPREL(TIME,  KHAT,  THAT,  TC) 

GIVEN  THE  MLE  OF  THE  COMPONENT  RELIABILITY  AND  THE  SAMPLE  SIZE, 

ENTER  THE  2-DIMENSIONAL  ARRAY  BI  AND  FIND  THE  BIAS  OF  THE 
ESTIMATOR 

DO  50  I  -  1,  NSMSZl 

LSMSZ  -  I 

IF  (RNSAM  .LE.  SMSZ(I-»-l))  GO  TO  60 
50  CONTINUE 

GO  TO  90 

60  IF  (RCHAT(J)  .LT.  .5)  GO  TO  100 

DO  70  I  -  1,  NRLBSl 

LRLBS  •»  I 

IF  (RCHAT(J)  .LE.  RLBS(I+1))  GO  TO  80 
70  CONTINUE 

GO  TO  90 

80  CALL  IBCICU  (BI,  ISD,  RLBS,  NRLBS,  SMSZ,  NSMSZ,  LRLBS,  LSMSZ, 

1  C,  WK,  lER) 

CALL  IBCEVU  (RLBS,  NRLBS,  SMSZ,  NSMSZ,  LRLBS,  LSMSZ,  C, 

1  RCHAT(J),  RNSAM,  BIAS,  lER) 

IF  (LRLBS  .LT.  NRLBS)  GO  TO  110 
90  BIAS  -  0. 

GO  TO  110 

100  IF  (RNSAM  .EQ.  SMSZ(LSMSZ  +  1))  LSMSZ  -  LSMSZ  +  1 

BIAS  -  Bid,  LSMSZ) 


50 


o  o  on  non 


110  RCHAT(J)  -  RCHAT(J)  -  BIAS 

ARCHAT(NCOUNT,  J)  -  RCHAT(J) 

C 

C  GIVEN  THE  UNBIASED  ESTIMATOR  OF  THE  COMPONENT  RELIABILITY, 

C  ENTER  THE  2-DIMENSIONAL  ARRAY  SIG  AND  FIND  THE  STANDARD 
C  DEVIATION  OF  THE  UNBIASED  ESTIMATOR 
Z  -  RCHAT(J) 

RCRLB  -  SQRT(Z  **  2  *  (ALOG(Z))  **  2  *  (1.109  -  .514  * 

1  ALOG(  -  ALOG(Z))  +  .608  *  (AL0G(  -  ALOG(Z)))  **  2)  /  RNSAM) 

DO  120  I  -  1,  NSMSZl 


LSMSZ 

»  I 

IF  (RNSAM 

.LE.  SMSZ(I+1)) 

GO  I 

120 

CONTINUE 

GO  TO  160 

130 

IF  (RCHAT(J) 

.LT.  .5)  GO  TO 

170 

130 

IF  (RCHAT (J) 

.LT.  .5)  GO  TO 

170 

DO  140  I 

-  1,  NRLBSl 

LRLBS 

-  I 

IF  (RCHAT(J)  .LE.  RLBS(I+1))  GO  TO  150 
140  CONTINUE 

GO  TO  160 

150  CALL  IBCICU  (SIG,  ISD,  RLBS,  NRLBS,  SMSZ,  NSMSZ,  LRLBS, 

1  LSMSZ,  C,  WK,  lER) 

CALL  IBCEVU  (RLBS,  NRLBS,  SMSZ,  NSMSZ,  LRLBS,  LSMSZ,  C, 

1  RCHAT(J),  RNSAM,  SIGMA,  lER) 

IF  (SIGMA  .GE.  RCRLB)  GO  TO  180 
160  SIGMA  -  RCRLB 

lER  -  0 

GO  TO  180 

170  IF  (RNSAM  .EQ.  SMSZ  (LSMSZ  +  D)  LSMSZ  -  LSMSZ  +  1 

SIGMA  -  SIG(1,  LSMSZ) 

IF  (SI(MA  .GE.  RCRLB)  GO  TO  180 
SIGMA  -  RCRLB 

180  CONTINUE 

FORM  A  VECTOR  OF  NDEV  SAMPLE  RELIABILITIES  WITH  MEAN  RCHAT 
AND  A  STANDARD  DEVIATION  OF  SIGMA 
NDEV  -  2000 

CALL  GGNML  (DSEED,  NDEV,  DEV) 

DO  190  I  -  1,  2000 

RC(I,  J)  -  RCHAT(J)  +  DEV(I)  *  SIGMA 
TRUNCATE  THE  NORMAL  DISTRIBUTION  IF  COMPONENT  RELIABILITY 
IS  GREATER  THAN  1 

IF  (RC(I,J)  .LE.  1.)  GO  TO  190 
RC(I,  J)  -  1. 

NTRUNC  -  NTRUNC  +  1 
190  CONTINUE 

GO  TO  290 


IF  FAILS  TO  CONVERGE,  START  WITH  ANOTHER  KHAT  AND  TRY  AGAIN 


200 

PRINT 

449,  NCOUNT,  J,  NSTART,  KHAT,  THAT,  ITMAX, 

lER 

NSTART 

-  NSTART  +  1 

GO  TO 

(  210,  220,  230,  240,  250,  260,  270,  280), 

NSTART 

210 

START 

••  1. 

51 


220 

GO  TO 
START 

30 

■ 

.5 

230 

GO  TO 
START 

30 

.9 

240 

GO  TO 
START 

30 

4.5 

250 

GO  TO 
START 

30 

a 

1.1 

260 

GO  TO 
START 

30 

1.5 

270 

GO  TO 
START 

30 

m 

3.5 

GO  TO 

30 

C  IF  MLES  FAIL  TO  CONVERGE  IN  8  ATTEMPTS,  PRINT  A 
C  MESSAGE  AND  GO  ON  TO  THE  NEXT  COMPONENT 
280  PRINT  459 

C 

290  CONTINUE 

Q  ***A*A*4c*A**A*4t**4tVk**4t4t**4t4r*A**4t*4t*A*A**A**************************** 

C 

C  FOR  EACH  OF  THE  4  SYSTEMS,  THE  DIFFERENT  COMPONENTS  ARE  COMBINED 
C  TO  YIELD  NDEV  SAMPLES  OF  THE  SYSTEM  RELIABILITY.  THESE  SAMPLES 
C  ARE  SEQUENCED  IN  ASCENDING  ORDER  AND  THEN  COUNTERS  KEEP  TRACK 
C  OF  WHEN  THE  99,  95,  90,  80,  70,  60,  AND  50  PERCENT  CONFIDENCE 
C  INTERVALS  CONTAIN  THE  TRUE  SYSTEM  RELIABILITY. 

C 

C  SYSTEM  1 

NMC  -  1 

IDIM  -  1 

CALL  RELl  (NMC,  TRC,  TRS,  IDIM) 

ATRS(l)  -  TRS(l) 

NMC  -  2000 

IDIM  -  2000 

CALL  RELl  (NMC,  RC,  RS,  IDIM) 

CALL  VSRTA  (RS,  IDIM) 

IF  (RS(20)  .LE.  TRS(l))  NS1C99  -  NS1C99  +  1 
IF  (RS(IOO)  .LE.  TRS(l))  NS1C95  -  NS1C95  +  1 

IF  (RS(200)  .LE.  TRS(l))  NS1C90  -  NS1C90  +  1 

IF  (RS(400)  .LE.  TRS(l))  NS1C80  -  NS1C80  +  1 

IF  (RS(600)  .LE.  TRS(l))  NSIC70  -  NS1C70  +  1 

IF  (RS(800)  .LE.  TRS(l))  NS1C60  -  NSIC60  +  1 

IF  (RS(IOOO)  .IE.  TRS(l))  NS1C50  -  NS1C50  +  1 
C 

C  SYSTEM  2 

NMC  -  1 

IDIM  -  1 

CALL  REL2  (NMC,  TRC,  TRS,  IDIM) 

ATRS(2)  -  TRS(l) 

NMC  -  2000 

IDIM  -  2000 

CALL  REL2  (NMC,  RC,  RS,  IDIM) 

CALL  VSRTA  (RS,  IDIM) 

IF  (RS(20)  .LE.  TRS(l))  NS2C99  -  NS2C99  +  1 
IF  (RS(IOO)  .LE.  TRS(l))  NS2C95  -  NS2C95  +  1 


52 


onon  no  no 


IF 

IF 

IF 

IF 

IF 


(RS(200)  .LE.  TRS(l))  NS2C90  -  NS2C90  +  1 

(RS(400)  .LE.  TRS(l))  NS2C80  -  NS2C80  +  1 

(RS(600)  .LE.  TRS(l))  NS2C70  -  NS2C70  +  1 

(RS(800)  .LE.  TRS(l))  NS2C60  -  NS2C60  +  1 

(RS(IOOO)  .LE.  TRS(l))  NS2C50  -  NS2C50  +  1 


SYSTEM  3 

NMC  -  1 

IDIM  -  1 

CALL  REL3  (NMC,  TRC,  TRS,  IDIM) 

ATRS(3)  -  TRS{1) 

NMC  -  2000 

IDIM  -  2000 

CALL  REL3  (NMC,  RC,  RS,  IDIM) 

CALL  VSRTA  (RS,  IDIM) 

IF  (RS(20)  .LE.  TRS(l))  NS3C99  -  NS3C99  +  1 
IF  (RS(IOO)  .LE.  TRS(l))  NS3C95  -  NS3C95  +  1 

IF  (RS(200)  .LE.  TRS(l))  NS3C90  -  NS3C90  +  1 

IF  (RS(400)  .LE.  TRS(l))  NS3C80  -  NS3C80  +  1 

IF  (RS(600)  .LE.  TRS(l))  NS3C70  -  NS3C70  +  1 

IF  (RS(800)  .LE.  TRS(l))  NS3C60  -  NS3C60  +  1 

IF  (RS(IOOO)  .LE.  TRS(l))  NS3C50  -  NS3C50  +  1 


SYSTEM  4 

NMC  -  1 

IDIM  -  1 

CALL  REL4  (NMC,  TRC,  TRS,  IDIM) 

ATRS(4)  -  TRS(l) 

NMC  =>  2000 

IDIM  =  2000 

CALL  REL4  (NMC,  RC,  RS,  IDIM) 

CALL  VSRTA  (RS,  IDIM) 

IF  (RS(20)  .LE.  TRS(l))  NS4C99  -  NS4C99  +  1 
IF  (RS(IOO)  .LE.  TRS(l))  NS4C95  -  NS4C95  +  1 

IF  (RS(200)  .LE.  TRS(l))  NS4C90  -  NS4C90  +  1 

IF  (RS(400)  .LE.  TRS(l))  NS4C80  -  NS4C80  +  I 

IF  (RS(600)  .LE.  TRS(l))  NS4C70  -  NS4C70  +  1 

IF  (RS(8O0)  .LE.  TRS(l))  NS4C60  -  NS4C60  +  1 

IF  (RS(IOOO)  .LE.  TRS(l))  NS4C50  -  NS4C50  +  1 
300  CONTINUE 


******************************************************************** 

******************************************************************** 


NPT  -  10000  *  NOLMC 

PRINT  469,  NTRDNC,  NPT 
RNOLMC  -  NOLMC 

C  FOR  SYSTEM  1,  DETERMINE  THE  99,  95,  90,  80,  70,  60,  AND  50 

C  PERCENT  CONFIDENCE  LIMIT  COVERAGE  OF  THE  TRUE  SYSTEM  RELIABILITY 
PS1C99  -  NS1C99  /  RNOLMC 

PS1C95  -  NS1C95  /  RNOLMC 

PS1C90  -  NS1C90  /  RNOLMC 

PS1C80  -  NS1C80  /  RNOLMC 

PS1C70  -  NS1C70  /  RNOLMC 


i 
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PS1C60  -  NS1C60  /  RNOLMC 
PS1C50  -  NS1C50  /  RNOLMC 

C  FOR  SYSTEM  2,  DETERMINE  THE  99,  95,  90,  80,  70,  60,  AND  50 
C  PERCENT  CONFIDENCE  LIMIT  COVERAGE  OF  THE  TRUE  SYSTEM  RELIABILITY 


PS2C99 

m 

NS2C99 

/ 

RNOLMC 

PS2C95 

m 

NS2C95 

/ 

RNOLMC 

PS2C90 

m 

NS2C90 

/ 

RNOLMC 

PS2C80 

m 

NS2C80 

/ 

RNOLMC 

PS2C70 

m 

NS2C70 

/ 

RNOLMC 

PS2C60 

m 

NS2C60 

/ 

RNOLMC 

PS2C50 

m 

NS2C50  / 

RNOLMC 

C  FOR  SYSTEM  3,  DETERMINE  THE  99,  95,  90,  80,  70,  60,  AND  50 
C  PERCENT  CONFIDENCE  LIMIT  COVERAGE  OF  THE  TRUE  SYSTEM  RELIABILITY 


PS3C99 

m 

NS3C99 

/ 

RNOLMC 

PS3C95 

m 

NS3C95 

/ 

RNOLMC 

PS3C90 

m 

NS3C90 

/ 

RNOLMC 

PS3C80 

m 

NS3C80  / 

RNOLMC 

PS3C70 

m 

NS3C70  / 

RNOLMC 

PS3C60 

m 

NS3C60 

/ 

RNOLMC 

PS3C50 

m 

NS3C50 

/ 

RNOLMC 

C  FOR  SYSTEM  4,  DETERMINE  THE  99,  95,  90,  80,  70,  60,  AND  50 
C  PERCENT  CONFIDENCE  LIMIT  COVERAGE  OF  THE  TRUE  SYSTEM  RELIABILITY 
PS4C99  -  NS4C99  /  RNOLMC 

PS4C95  -  NS4C95  /  RNOLMC 

PS4C90  -  NS4C90  /  RNOLMC 

PS4C80  -  NS4C80  /  RNOLMC 

PS4C70  -  NS4C70  /  RNOLMC 

PS4C60  -  NS4C60  /  RNOLMC 

PS4C50  -  NS4C50  /  RNOLMC 


PRINT  479 

PRINT  489,  ATRS 

(1). 

PS1C99, 

PS1C95, 

PS1C90, 

PS1C80, 

PS1C70 

PS1C60,  PS1C50 
PRINT  499 

PRINT  489,  ATRS 

(2), 

PS2C99, 

PS2C95, 

PS2C90, 

PS2C80, 

PS2C70 

PS2C60,  PS2C50 
PRINT  509 

PRINT  489,  ATRS 

(3), 

PS3C99, 

PS3C95, 

PS3C90, 

PS3C80, 

PS3C70 

PS3C60,  PS3C50 
PRINT  519 

PRINT  489,  ATRS 

(4), 

PS4C99, 

PS4C95, 

PS4C90, 

PS4C80, 

PS4C70 

PS4C60,  PS4C50 

C  *4t*4tAA4t*AA4t4tAA*A*4i*4tA*4t***4t*A****4tA4i*********AA****ilt*A*****4t**4iAik*4t* 
Q  ****4tA4l*4t4t*A*4t4i*4tA**A4t*A4l4t*A4tA*A**A*****A*A***4tA*A4tAAAA****4t*A)lrA**** 

C 

C  THIS  PROGRAM  SEOIENT  RUNS  A  KOLMOGOROV-SMIRNOV  TEST  OF  THE 
C  ASSUMPTION  THAT  THE  COMPONENT  RHAIS  ARE  NORMALLY  DISTRIBUTED 
C 

DO  360  J  -  1,  5 

Z  -  TRC(J) 

RCRLB  -  SQRT(Z  **  2  *  (ALOG(Z))  **  2  *  (1.109  -  .514  * 

1  ALOG(  -  ALOG(Z))  +  .608  *  (ALOG(  -  ALOG(Z)))  **  2)  /  RNSAM) 

DO  310  I  -  1,  NSMSZl 

LSMSZ  -  I 
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IF  (RNSAM 

.LE.  SMSZ(I+1)) 

GO  TO  320 

310 

CONTINUE 

GO  TO  350 

320 

DO  330  I 

-  1,  NRLBSl 

LRLBS 

-  I 

IF  (TRC(J) 

.LE.  RLBS(I+1)) 

GO  TO  340 

330 

CONTINUE 

GO  TO  350 

340 

CALL  IBCICU 

(SIG,  ISD,  RLBS, 

NRLBS,  SMSZ,  NSMSZ,  LRLBS,  LSMSZ, 

1  C,  WK,  lER) 


CALL  IBCEVU  (RLBS,  tniLBS,  SMS2,  NSMSZ,  LRLBS,  LSMSZ,  C,  TRC(J), 
1  RNSAM,  SIGMA,  lER) 

IF  (SIGMA  .GE.  RCRLB)  GO  TO  360 
350  SIQIA  -  RCRLB 
lER  -  0 

PRINT  *,  "Z-",  Z,  "  SIGMA-",  SIGMA 
360  RCHSIG(J)  -  SIGMA 

PRINT  529,  (TRC(J),  J-1,  5) 

PRINT  529,  (RCHSIG(J),  J-1,  5) 

PRINT  *,  "  " 

ALPHA  -  .05 

TALPHA  -  1.36  /  SQRT(RNOLMC) 

IF  (NOLMC  .LE.  40)  TALPHA  -  QNTKS  (NOLMC) 

CALL  TNORKS  (ARCHAT,  NOLMC,  TRC,  RCHSIG,  ALPHA,  TALPHA) 

C 

C  ******************************************************************** 
Q  ******************************************************************** 

c 

STOP  "FORAWHILE" 


409  FORMAT  (  IHI  ) 

419  FORMAT  (  "  ***************************************************" 


1 

"**********<•  ^ 

/  .  " 

*",  T62, 

/  . 

"  *",  T62, 

2 

ft 

f 

T25, 

"SAMPLE 

"  "SIZE  -  ", 

13, 

T62,  /  , 

3 

T62, 

/ 

If 

9  9 

T22,  "MONTE 

CARLO 

"  "SIZE  -  ",  13 

4 

T62, 

/ 

II  All 

9  ^  9 

T62,  "*",  / 

tl 

9 

T62, 

5  "  ************************************************************" 


6 

'»*" 

f 

/////// 

) 

429 

FORMAT 

(  5 (IX,  "COMPONENT  " 

9 

11,  /  ,  6X, 

"K 

-  ",  F4.2,  / 

9 

1 

6X, 

"THETA  -  ",  F5.0,  /  , 

6X,  "C  -  ",  F2 

:.o. 

/  .  6X, 

2 

"RELIABILITY  -",  F7.5,  / 

/ 

)  ) 

439 

FORMAT 

(  15(1X,  F7.5)  ) 

449 

FORMAT 

{  !  ,  "  **********" 

9 

/  ,  "  MC-", 

14, 

4X,  "J  -", 

12, 

1 

4X, 

"NSTART  -",12,  /  , 

9f 

KHAT-"  E13.6, 

5X 

"THAT-"  E13.6 

/ 

2 

"  ITERATIONS-"  13,  5X  "lER- 

ft 

13  /  "  **********"  ) 

459 

FORMAT 

(  /  '■  *******************"  II" 

DID 

NOT  CONVERGE 

IN" 

1  "8  ATTEMPTS  WITH  DIFFERENT  STARTING  KHATS"  /  "  THEREFORE  GO" 

2  "ING  ON  TO  THE  NEXT  COMPONENT"  II"  ********************” 

3  ) 

469  FORMAT  (  /  /  ,  "  THERE  WERE",  17,  "  HIGH  TRUNCATIONS  OUT  OF" 

1  ,  18,  "  RELIABILITY  DEVIATES",  /  /  ) 

479  FORMAT  (  /  "  *****  SYSTEM  1  *****"  /  "  (3  COMPONENTS  IN  SERI" 

1  "ES)"  ) 
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r 


489  FORMAT  (  /  ,  "  TRUE  SYSTEM  RELIABILITY  F7.5,  / 

1  "  THE  99  PERCENT  CONFIDENCE  INTERVAL  COVERED  ",  F6.4, 

2  "  OF  THE  RUNS",  /  ,  "THE  95  PERCENT  CONFIDENCE  INTERVAL  COV" 

3  "ERED  ",  F6.4  "OFTHERUNS",  /  ,  "THE  90  PERCENT  CONFIDENCE  I" 

4  "NTERVAL  COVERED  "  F6.4,  "  OF  THE  RUNS",  /  ,  "THE  80  PERCE" 

5  "NT  CONFIDENCE  INTERVAL"  "  COVERED  "  F6.4  "OFTHERUNS”,  /  , 

6  "  THE  70  PERCENT  CONFIDENCE  "  "INTERVAL  COVERED  ",  F6.4, 

7  "  OF  THE  RUNS",  /  ,  "THE  60  PERCENT  "  "CONFIDENCE  INTERVAL" 

8  "  COVERED  ",  F6.4,  "  OF  THE  RUNS",  /  ,  "  THE",  "  50  PERCEN" 

9  "T  CONFIDENCE  INTERVAL  COVERED  ",  F6.4,  "OF  THE  RUNS",  /  / 

9  /  ) 


499 

FORMAT  ( 

/  " 

*****  SYSTEM 

2 

*****” 

/ 

"  (1  COMPONENT  IN  SERIE" 

1 

"S  WITH 

2  " 

"IN  PARALLEL)" 

) 

509 

FORMAT  ( 

/  " 

*****  SYSTEM 

3 

*****” 

/ 

"  (3  CCMPONENTS  IN  PARA" 

1 

"LLED" 

) 

519 

FORMAT  ( 

/  " 

*****  SYSTEM 

4 

*****" 

/ 

"  (A  5-COMPONENT  COMPLE" 

1 

"X  NETWORK)" 

) 

i 

529 

FORMAT  ( 

5(1X, 

F10.8,  5X)  ) 

END 

FUNCTION  CMPREL  (TIME,  K,  THETA,  C) 

REAL  K 

CMPREL  -  EXP (  -  ((TIME  -  C)  /  THETA)  **  K) 

RETURN 

END 
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FUNCTION  RLKLHD  (R,  RK,  T,  C) 

DIMENSION  R  (205,5) 

J  -  R(205,  1) 

NSAM  -  R(204,  J) 

RLKLHD  -  1 . 

DO  10  I  -  1,  NSAM 
X  -  RLKLHD 

10  RLKLHD  -  X  *  RK  *  (R(I,  J)  -  C)  **  (RK  -  1.)  /  T  **  RK  * 
1  EXP(  -  (((R(I,  J)  -  C)  /  T)  **  RK)) 

RETURN 

END 
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FUNCTION  F  (X,  K,  R) 

DIMENSION  R  (205,5),  X  (2) 

REAL  KHAT 

J  -  R(205,  1) 

NSAM  -  R(204,  J) 

KHAT  -  X(l) 

THAT  -  X(2) 

SUMXIK  -  0. 

SUMLXI  -  0. 

SUMFR  -  0. 

IF  (THAT  .LE.  0.  .OR.  KHAT  .GT.  15.)  RETURN 
DO  10  I  -  1,  NSAM 

SUMXIK  -  SUMXIK  +  R(I,  J)  **  KHAT 

10  SUMLXI  -  SUMLXI  +  ALOG(R(I,  J)) 

DO  20  I  -  1,  NSAM 

20  SUMFR  -  SUMFR  +  (R(I,  J)  /  THAT)  **  KHAT  *  ALOG(R(I,  J)  / 

1  THAT) 

IF  (K  .EQ.  2)  GO  TO  30 

F  -  THAT  -  (SUMXIK  /  NSAM)  **  (1.  /  KHAT) 

RETURN 

30  F  »  NSAM  /  KHAT  -  NSAM  *  ALOG(THAT)  +  SUMLXI  -  SUMFR 

RETURN 
END 
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SUBROUTINE  RELl  (NMC,  R,  RS,  IDIM) 

C  RELl  DETERMINES  THE  SYSTEM  RELIABILITY  OF  3  COMPONENTS  IN  SERIES 
DIMENSION  R  (IDIM, 5),  RS  (IDIM) 

DO  10  I  »  1,  NMC 

10  RS(I)  -  R(I,  1)  *  R(I,  2)  *  R(I,  3) 

RETURN 

END 
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SUBROUTINE  REL2  (NMC,  R,  RS,  IDIM) 

C  REL2  DETERMINES  THE  SYSTEM  RELIABILITY  OF  1  COMPONENT  IN 
C  SERIES  WITH  2  IN  PARALLEL 

DIMENSION  R  (IDIM, 5),  RS  (IDIM) 

DO  10  I  -  1,  NMC 

10  RS(I)  -  R(I,  1)  *  (1.  -  (1.  -  R(I,  2))  *  (1.  -  R(l,  3))) 
RETURN 
END 
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SUBROUTINE  REL3  (NMC,  R,  RS,  IDIM) 

C  REL3  DETERMINES  THE  SYSTEM  RELIABILITY  OF  3  COMPONENTS  IN  PARALLEL 
DIMENSION  R  (IDIM, 5),  RS  (IDIM) 

DO  10  I  -  1,  NMC 

10  RS(I)  -  1.  -  (1.  -  R(I,  D)  *  (1.  -  R(I,  2))  *  (1.  -  R(I, 
1  3)) 

RETURN 

END 
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SUBROUTINE  REL4  (NMC,  R,  RS,  IDIM) 

C  REL4  DETERMINES  THE  SYSTEM  RELIABILITY  OF  A  5  COMPONENT 
C  COMPLEX  NETWORK 

DIMENSION  R  (IDIM, 5),  RS  (IDIM) 

DO  10  I  -  1,  NMC 

RS(I)  -  R(I,  1)  *  (1  -  (1  -  R(I,  2))  *  (1  -  R(I,  5)  *  (1  - 
1  (1  -  Rd.  3))  *  (1  -  R(I,  4))))) 

10  CONTINUE 
RETURN 
END 


V 

'  1 


J 
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SUBROUTINE  TNORKS  (ARCHAT,  NOLMC,  TRC,  RCHSIG,  ALPHA,  TALPHA) 
DIMENSION  ARCHAT  (2000,5),  CHMX  (5),  RCHSIG  (5),  TRC  (5) 

RNOLMC  -  NOLMC 

C  ORDER  EACH  COLUMN  (A  VECTOR  CONTAINING  THE  COMPONENT  RHATS) 

C  FROM  THE  SMALLEST  RHAT  TO  THE  LARGEST 
DO  10  J  -  1,  5 

10  CALL  VSRTA  (ARCHAT(1,  J) ,  NOLMC) 

C 

C  **********  for  each  COMPONENT  ********** 

DO  40  J  -  1,  5 
C 

C  FIND  THE  DEVIATIONS  OF  THE  EMPERICAL  DISTRIBUTION  FROM  THE 
C  CUMULATIVE  PROBABILITY  DISTRIBUTION  OF  THE  NORMAL 
AVRCH  -  0. 

DO  20  I  -  1,  NOLMC 

C  AVRCH  SUMS  THE  RHATS  FOR  EACH  COMPONENT.  LATER  IT  WILL  BE 
C  DIVIDED  BY  THE  MONTE  CARLO  SIZE  TO  GET  AN  AVERAGE  RHAT. 

AVRCH  -  AVRCH  +  ARCHAT (I,  J) 

SD  -  (ARCHATd,  J)  -  TRC(J))  /  RCHSIG(J) 

CALL  MDNOR  (SD,  X) 

Z  «  I 

CHLO  -  ABS((Z  -  1.)  /  RNOUIC  -  X) 

CHHI  »  ABS(Z  /  RNOLMC  -  X) 

IF  (I  .EQ.  1)  CHMX  (J)  -  CHLO 
IF  (CHLO  .GT.  CHMX(J))  CHMX  (J)  »  CHLO 
IF  (CHHI  .GT.  CHMX(J))  CHMX  (J)  »  CHHI 
PRINT  *,  J,  I,  ARCHAT  (I,  J) ,  X,  CHLO,  CHHI 
20  CONTINUE 

C  PRINT  THE  RESULTS  OF  THE  K-S  TEST  OF  NORMALITY 
AVRCH  -  AVRCH  /  RNOLMC 

BIAS  -  AVRCH  -  TRC(J) 

ESD  -  0. 

DO  30  I  -  1,  NOLMC 

DIFF  -  ARCHATd,  J)  -  AVRCH 

ESD  «■  ESD  +  DIFF  *  DIFF 

30  CONTINUE 

ESD  -  SQRT(ESD  /  (RNOLMC  -  1.)) 

PRINT  109,  J,  TRC  ( J) ,  RCHSIG  ( J) ,  AVRCH,  ESD,  BIAS,  CHMX  ( J) , 
1  TALPHA 

IF  (CHMX(J)  .GT.  TALPHA)  PRINT  119, ALPHA 
IF  (CHMX(J)  .LE.  TALPHA)  PRINT  129, ALPHA 
40  CONTINUE 

Q  **************************************** 

C 

RETURN 

C 

C 

109  FORMAT  (  /  ,  IX,  "K-S  TEST  OF  NORMALITY  FOR  RHATS  OF  COMPONEN" 

1  "T  ",  II,  /  ,  5X,  "TRUE  COMPONENT  RELIABILITY  -",  F10.8,  /  , 

2  5X,  "SIGMA  -",  F10.8,  /  ,  5X,  "AVERAGE  RHAT  F10.8,  /  , 

3  5X,  "EMPERICAL  STANDARD  DEVIATION  F10.8,  /  ,  5X, 

4  "BIAS  -",  F10.8,  /  ,  5X,  "TEST  STATISTIC  -",  F5.3,  8X, 

5  "TALPHA  -",  F5.3  ) 

119  FORMAT  (  5X,  "REJECT  AT  ALPHA  F4.2  ) 

64 


129  FORMAT  (  5X 


"FAIL  TO  REJECT  AT  ALPHA 


,11 


F4.2  ) 


C 


END 


I 
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